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coherence  measures  can  then  be  used  to  test  the  effectiveness  of  such  factorization. 
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19. 


Following  a  brief  discussion  of  spectral  factorabilily  and  motivations  for  broadband  analysis,  the  report 
is  subdivided  into  four  main  sections.  In  Section  1.0,  we  examine  the  spectral  factorability  of  Pn  and  Lg  phases 
recorded  at  NORESS  and  the  3C  subarray  of  NORSAR.  We  find  that  closely  located  events  can  be  subdivided  into 
distinct  factorable  groups  with  coherence  measures  indicating  high  effectiveness  of  factorization.  Events  not 
within  such  groups  do  not  factor  efficiently,  and  thus,  result  in  a  low  coherence  measure.  The  factorizable  groups 
are  probably  events  which  are  close  in  location,  depth  and  source  mechanism.  Within  factorizable  groups,  fine 
details  of  regional  waveforms  can  be  reconstructed  from  source  and  site  factors.  This  is  the  first  utilization  of  the 
information  in  the  details  of  regional  waveforms  we  know  of;  even  the  most  ambitious  direct  modeling  efforts  can 
reproduce  only  envelope  shapes  and  a  few  main  features,  such  as  relative  phase  amplitudes  of  regional 
seismograms.  Factorization  studies  can  thus  result  in  better  location  and  the  detection  of  anomalous  events  with 
differing  source  mechanisms  or  depths  (such  as  decoupled  large  explosions  in  some  evasion  scenarios).  Methods 
are  described  for  diagnosing  regional  waveform  differences  due  to  relative  location  in  source  time  functions, 
source  depths  and  mechanisms  using  standard  linear  system  concepts.  Neither  the  source  time  functions  nor  the 
Green’s  functions  need  be  known  for  such  analyses. 

Section  2.0  re-examines  some  aspects  of  estimating  secondary  arrivals  (pP,  spall,  etc.)  in  teleseismic  P 
waves  and  applying  corrections  to  m^.  A  popular  model  (P+pP  model)  often  used  for  interpreting  P  waves  from 
nuclear  explosions  and  deriving  pP  parameters  consists  of  a  site  function  convolved  with  a  P  and  pP  pulse 
sequence  and  some  explosion  source  pulse  (Lay,  1985;  Murphy,  1989).  Frequency  domain  analyses  of  data 
variance  show  that  this  model  cannot  explain  a  significant  portion  of  the  energy  (30  -  40%)  in  t.he  P  waves  from 
Pahute  Mesa  events  for  which  such  analysis  methods  were  previously  applied.  Spectral  factorization,  which  does 
not  assume  the  P+pP  model,  on  the  other  hand,  accounts  for  90  -  95%  of  the  energy  in  both  Pahute  Mesa  and 
Kazakh  events.  Therefore,  the  complex  source  waveforms  often  seen  in  multi-channel  deconvolved  P  waves  from 
nuclear  explosions  at  both  Pahute  Mesa  and  the  Kazakh  test  sites  are  required  for  explaining  the  data  and  reflect  the 
actual  complexity  of  P  waves  radiated  from  nuclear  explosions.  It  is  shown  that  magnitude-yield  relationships 
appropriate  to  the  generic  pP  characteristics,  derived  from  multi-chaitnel  deconvolutions  of  nuclear  explosions 
(and  appropriate  t*),  are  similar  to  those  directly  derived  by  regression  analyses  of  mb  data  and  published  yields  of 
nuclear  explosions  for  Pahute  and  Kazakh  test  sites  by  Jih  (1990). 

Section  3.0  explores  ways  to  perform  spatio-temporal  source  inversion  using  waveform  data.  The 
formalisms  used  are  generalizations  of  F-K  algorithms.  Shumway's  F-K  method,  compounded  over  various 
frequencies,  is  used  mostly  in  this  part  of  the  report.  Unlike  many  of  the  F-K  statistics  in  use,  the  statistics  for 
each  frequency  has  a  known,  non-central  F,  statistical  distribution.  Assuming  that  the  compounded  statistics  is 
approximately  normal,  confidence  bounds  on  the  various  parameters  derived  could  be  estimated.  An  alternative 
statistic  having  a  central  distribution  has  also  been  derived.  We  explored  the  resolution  of  inversion  of  long 
period  body  waveforms  and  concluded  that  it  is  extremely  poor  and  we  suggest  that  such  methodologies  must  be 
replaced  by  more  precise  methods.  The  concepts  are  demonstrated  by  analyses  of  actual  data  from  nuclear 
explosions  and  three  earthquakes;  the  1966  El  Golfo,  the  1976  Kazakh  and  the  1988  Armenian  events.  Much 
needs  to  be  done  for  further  refinement  of  the  methodology  we  have  developed. 

In  the  last  part  of  the  report  we  describe  an  approach  for  S/N  ratio  enhancement,  using  the  application  of 
site  waveform-equalization,  together  with  adaptive  filtering.  This  combines  the  reduction  of  the  beam  loss  with 
die  advantages  of  optimum-filtering.  The  method  was  applied  to  both  regional  and  teleseismic  array  data.  We 
lour.u  that,  typically,  we  were  able  to  obtain  a  4  -  6  dB  improvement  relative  to  conventional  beaming.  This  S/Tl 
r.itio  improvement  is  significant  enough  to  be  pursued  further,  despite  the  complexity  of  the  processing. 
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Figure  1:  Decomposition  and  reconstruction  of  a  matrix  of  seismograms,  using  site  and  source 
spectral  factor  representation. 

Figure  2;  Scenarios  for  which  spectral  factorability  may  be  valid,  a)  Teleseismic  P  from  a  limited 
source  region  recorded  at  a  seismic  array,  b)  Network  recordings  of  groups  of  events  with  similar 
radiation  patterns. 

Figure  3:  Flow  diagram  of  the  spectral  factorization  process  utilized  in  this  report. 

Figure  4:  Map  of  Scandinavia  with  the  location  of  NORESS  and  Titania,  Blasjo  and  the  Estonian 
mines  studied. 

Figure  5:  a)  Selected  reconstructed  and  data  traces  of  Pn  arrivals  from  the  Titania  1  event,  b)  The 
same  for  the  Titania4  event.  The  designation  C2  refers  to  the  second  sensor  in  the  C  ring  of 
NORESS,  etc. 

Figure  6:  a)  Selected  reconstructed  and  data  traces  of  Lg  arrivals  from  the  E4-1  event,  b)  The 
same  for  the  E4-4  event.  The  designation  C2  refers  to  the  second  sensor  in  the  C  ring  of 
NORESS,  etc. 

Figure  7:  Spectra  for  various  regional  arrivals  from  the  Titania  mine.  Note  that  the  spectral 
modulation  is  similar,  regardless  of  the  type  of  arrival. 

Figure  8:  Inter-event  coherences  for  Pn  arrivals  from  the  Titania  mine. 

Figure  9a:  Wiener  filtering  results  for  transforming  Event  Titania4  into  Event  Titanial. 

Figure  9b:  Wiener  filtering  results  for  transforming  Event  Titania4  into  Event  TitaniaS. 

Figure  10:  Cross-correlations  of  the  site -equalized  traces  of  Pn  between  sensors  D1  and  D5  for 
selected  Titania  events. 

Figure  11:  Site-equalized  Pn  traces  for  selected  Titania  events  between  sensors  D1  and  D5. 

Figure  12:  Cross-correlation  of  the  site-equalized  traces  and  the  equalized  traces  of  Lg  between 
sensors  D1  and  D5  for  the  E9  event.  Note  the  excellent  similarity  between  the  cross-equalized 
traces. 

Figure  13:  Results  of  site-equalization  when  spectral  factors  from  joining  factorization  of  El  and 
E4  Events  were  utilized.  Although  the  cross-correlations  still  have  peaks  at  the  proper  alignment 
with  reduced  correlation  coefficients,  (a)  the  efficiency  of  site-equalization  is  greatly  degraded  as 
evidenced  by  the  poor  similarity  of  the  site-equalized  traces  (b). 

Figure  14:  Site-averaged  coherence  results  from  the  joint  factorization  of  three  E9  events  and  one 
E4  event. 

Figure  15:  Site-averaged  coherence  results  from  the  joint  factorization  of  four  E9  events  with 
reconstructions  of  one  event  being  much  less  similar  to  the  original  data  than  the  others. 

Figure  16:  Envelopes  of  band-pass  filtered  traces  Type  1  of  E9  Event. 
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Figure  17;  Envelopes  of  band-pass  filtered  traces  Type  II  of  Event  E9.  Even  though  these  events 
seem  to  be  noisier  than  those  if  Figure  16,  the  S/N  ratio  for  for  Lg  for  these  events  is  about  the 
same  as  for  E9. 

Figure  18:  Coherences  demonstrating  the  high  factorization  efficiency  -  90  to  95%  of  the  energy 
accounted  for,  across  the  3  C  subarray  of  NORSAR,  of  the  Pn  phases  from  a  group  of  closely 
located  FENNOLORA  profiling  shots.  What  is  noteworthy  about  this  is  the  relatively  large 
spacing  about  5  km,  of  the  subarray  sensors. 

Figure  19:  Methods  for  testing  the  efficiency  (portion  of  energy  accounted  for),  a)  By  spectral 
factorization,  b)  By  the  P+pP  model. 

Figure  20:  Some  typical  deconvolution  results  for  various  test  sites. 

Figure  21 :  Data  and  reconstructions  using  the  spectral  factorization  method. 

Figure  22:  Waveform  inter-correlation  results  at  NORSAR  using  Murphy's  (1989)  pP  parameters 
for  Pahute  Mesa  events.  The  waveforms  at  the  same  sensors  should  be  identical,  but  they  are  not 
Indicating  that  the  P-i-pP  model  and/or  the  pP  parameters  are  invalid. 

Figure  23:  Portion  of  P  wave  energy  in  Pahute  Mesa  events  accounted  for  at  NORSAR  by  Lay’s 
(top  row),  Murphy’s  pP  parameters  (second  row)  all  shown  by  dashed  lines.  The  energy 
accounted  for  in  Pahute  Mesa  events  by  spectral  factorization  (multi-channel  deconvolution)  is 
shown  by  solid  lines  in  the  first  two  rows.  The  bottom  graph  shows  the  average  value  of  the  same 
for  Kazakh  events  obtained  by  spectral  factorization. 

Figure  24:  Simulated  WWSSN  mb  for  the  three  test  sites  obtained  by  using  the  source  model  of 
von  Seggem  and  Blandford  (1972),  t*  derived  from  P  spectra  and  pP  parameters  derived  from 
deconvolution.  Note  the  large  NTS-Kazakh  bias  and  the  Degelen  negative  offset  relative  to 
Shagan.  The  Ms-mb  lines  have  a  pronounced  curvature  and  the  biases  increase  with  decreasing 
yiela.  The  two  sets  off  curves  are  for  stations  located  in  shield  and  tectonic  areas  respectively. 
Identifical  source  media  were  assumed  for  simplicity. 

Figure  25:  Simulated  WWSSN  mb  for  the  three  test  sites  obtained  by  using  the  source  model  of 
Mueller  and  Murphy  (1971),  t*  derived  from  P  spectra  and  pP  parameters  derived  from 
deconvolution.  Note  the  large  NTS-Kazakh  bias  and  the  Degelen  Negative  offset  relative  to 
Shagan.  The  Ms-mb  lines  have  little  curvature  and  the  biases  increase  only  slightly  with 
decreasing  yield.  The  latter  features  are  more  in  agreement  with  the  data.  Identical  source  media 
were  assumed  for  simplicity. 

Figure  26:  Jih’s  regression  results  for  NTS  and  Kazakh  (Degelen  and  Shagan)  events  using 
maximum  likelihood  analyses  of  actual  and  simulated  P  waves  (recorded  at  other  stations 
transformed  to  WWSSN  response).  Note  the  offset  of  Degelen  relative  to  Shagan  and  the  large 
NTS- Shagan  difference  that  increase  with  decreasing  mb-  Note  the  similarities  of  these  hnes  to  the 
general  features  in  Figure  B1  and  especially  to  those  in  Figure  25. 

Figure  27:  Waveforms  used  for  computing  the  data  synthetic  coherences  in  Figure  3.  The 
different  traces  are  shown  be'  iw  each  pair.  The  two  numbers  next  to  each  pair  are  the  reference 
number  in  Table  II  and  the  time  domain  correlation  coefficient. 

Figure  28:  Averaged  spectra  of  a  sample  of  long  period  data  and  synthetic  seismogram  pairs 
randomly  selected  from  the  literature. 
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Figure  29:  Ensemble-averaged  coherences  of  our  sample  of  data  synthetic  pairs  defining  the 
“average  fit”.  The  95%  limit  shows  the  level  below  which  the  coherences  are  not  significantly 
different  from  zero. 

Figure  30:  Perspective  views  of  the  variations  of  the  frequency-compounded  F  statistics  over  the 
frequency  range  of  0.02-0.15  Hz  for  a  pair  of  double-couple  point  sources  with  equal  power. 
Uniform  random  noise  was  added  to  simulate  a  fit  with  an  average  correlation  coefficient  of  0.9. 
The  peak  value  occurs  between  the  two  sources  and  the  two  sources  were  not  resolved.  The  F 
statistics  were  computed  for  various  assumed  source  depths  and  the  sides  of  the  source  region 
views  are  250  km  long.  The  maximum  value  of  the  statistic  and  its  location  in  gridpoint 
coordinates  from  the  left  lower  comer  (5  km  spacing  was  used  in  both  x  and  y)  are  given  beside 
each  view.  The  gridpoint  x-y  coordinates  of  the  two  syntlietic  sources  were  at  27,  27  and  31,31 
and  the  depths  were  at  3  and  9  km  respectively. 

Figure  31:  Perspective  views  of  the  variations  of  the  frequency-compounded  F  statistics  over  the 
frequency  range  of  0.02-0.35  Hz  for  a  pair  of  double-couple  point  sources  the  F  values  were 
computed  for  various  assumed  source  depths.  Random  noise  with  16%  of  the  signal  was  added 
throughout  the  band.  The  peak  values  occur  nearer  to  the  two  sources  and  the  two  sources  were 
marginally  resolved  (for  ad^tional  details  of  this  display  see  the  caption  of  Figure  4). 

Figure  32:  Perspective  views  of  the  variations  of  the  frequency-compounded  F  statistics  over  the 
frequency  range  of  0.15-0.35  Hz,  thus  deleting  the  low  frequencies,  for  a  pair  of  double-couple 
point  sources  the  F  values  were  computed  for  various  assumed  source  depths.  Random  noise  with 
16%  of  the  signal  was  added  throughout  the  band.  The  peak  values  occur  near  to  the  two  sources 
and  the  two  sources  were  well  resolved  (for  additional  details  of  this  display  see  the  caption  of 
Figure  4). 

Figure  33:  Perspective  views  of  the  variations  of  the  positive  parts  of  the  simulated  site-averaged 
time  domain  data  vs.  probe  vector  (data  synthetic)  correlation  coefficient  over  the  frequency  range 
of  0.02-0.15  Hz  for  a  pair  of  double-couple  point  sources  the  F  values  were  computed  for  various 
assumed  source  depths.  The  negative  parts  of  the  correlation  coefficient  were  set  to  zero 
arbitrarily,  since  only  peak  positive  values  are  of  interest.  Random  noise  was  added  throughout 
the  band  to  give  a  peak  value  near  0.9.  The  two  sources  were  rot  resolved  (for  additional  details 
of  this  display  see  the  caption  of  Figure  4). 

Figure  34:  a)  F  statistics  computed  for  the  set  of  synthetic  seismograms  appropriate  to  the  El 
Golfo  crustal  and  point  source  model  at  10  km  depth  as  fitted  to  a  wide  range  of  depths.  Five 
percent  random  white  noise  was  added  to  avoid  singularity  in  the  solution,  llie  dashed  lines  are 
approximate  95%  confidence  limits  to  the  F  values.  The  maximum  unambiguously  defines  a  depth 
of  10  km.  b)  F  statistics  computed  using  the  actual  P  waves  from  the  El  Golfo  earthquake  and 
fitted  to  a  wide  range  of  depths.  The  dashed  lines  are  approximate  95%  confidence  limits.  The 
arrows  indicate  the  large  acceptable  depth  range  which  is  consistent  with  the  F  values. 

Figure  35:  Selected  data  and  synthetic  seismograms  appropriate  to  various  source  depths  for  El 
Golfo  earthquake.  The  “best”  time  domain  waveform  fits  occur  at  difference  depths  at  the  various 
stations. 

Figure  36a:  Matching  the  three  subevent  source  models  for  the  Armenian  earthquake  against 
themselves  in  the  0.1-1  i)  Hz  range  using  both  perspective  views  and  contour  plots.  Note  that  the 
distance  scales  are  different  on  the  horizontal  (EW)  and  the  vertical  (NS)  ax’ a)  Subevent  1. 
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Figure  36b:  Subevent  2. 


Figure  36c:  Subevent  3.  Perspective  views  are  on  to  and  contour  plots  of  the  same  surface  are 
below.  The  sizes  of  grids  are  127.5  km  x  127.5  km  with  51  subdivisions  in  each  direction. 

Figure  37a:  Matching  the  three  subevent  source  models  for  the  Armenian  earthquake  against  each 
other  in  the  0.1 -1.0  Hz  range,  a)  Model  2  against  Model  1. 

Figure  37b:  Model  3  against  Model  1. 

Figure  37c:  Model  3  against  Model  2.  Note  that  the  distance  scales  are  different  on  the  horizontal 
(EW)  and  the  vertical  (NS)  axis.  Perspective  views  are  on  top  and  contour  plots  of  the  same 
surface  are  below.  The  sizes  of  grids  are  127.5  km  x  127.5  Ion  with  51  subdivisions  in  each 
direction. 

Figure  38a:  Matching  the  three  subevent  source  models  for  the  Armenian  earthquake  against  P 
wave  data  form  eight  teleseismic  stations  in  the  0.1 -1.0  Hz  range.  Matching  with  a)  Subevent  1. 

Figure  38b:  Subevent  2. 

Figure  38c:  Subevent  3.  Note  that  the  dominant  feature  is  a  EW  trending  ridge.  Perspective 
views  are  on  top  and  contour  plots  of  the  same  surface  are  below.  The  size  of  grids  are  127.5  km 
X  127.5  km  with  51  subdivisions  in  each  direction. 

Figure  39:  Shumway’s  statistic  as  a  function  of  depth  computed  from  deconvolved  P  wave 
seismograms  of  the  December  25,  1975  Kazakh  nuclear  explosion. 

Figure  40:  Shumway’s  statistic  as  a  function  of  depth  computed  from  deconvolved  P  wave 
seismograms  of  the  June  29, 1977  Kazakh  nuclear  explosion. 

Figure  41:  Surface  wave  imaging  of  two  nuclear  explosions  at  the  Kazakh  test  site  (October  12, 
1980  on  top  and  September  14,  1980  on  the  bottom )  by  the  Capon  algorithm.  The  geographical 
coordinates  of  the  centers  of  the  print  plot  grids  are  50N-79E  and  the  events  are  plotted  relative  to 
that .  Filled  circles  denote  the  locations,  obtained  from  body  waves ,  as  listed  by  unclassified  files 
at  the  Center  for  Seismic  Studies,  triangles  show  maxima  of  the  Capon’s  statistics.  The  grids  are 
50  km  on  each  side,  and  the  tops  are  facing  North. 

Figure  42:  Deconvolved  waveforms  of  the  Kazakh  earthquake  at  three  UK  arrays.  Besides  the 
surface  reflections  pP  and  sP  there  is  an  early  phase  which  is  probably  a  Moho-converted  phase 
(Pooley  et  at  1983). 

Figure  43:  Depth  imaging  of  the  3/20/1076  Kazakh  earthquake. 

Figure  44:  Comparisons  of  noise  adaptation  with  out  (bottom)  and  with  (top)  the  site  response 
removed.  Comparing  the  shapes  and  levels  of  spectra  it  can  be  sent  hat  the  overall  noise  is 
significantly  reduced  by  the  adaptive  beaming  method  relatively  to  the  beams,  especially  below  2 
Hz.  The  adaptive  beams  quickly  adapt  to  the  noise,  and  after  a  few  oscillations  the  overall  noise 
levels  are  much  reduced  relative  to  the  beam  outputs. 

Figure  45:  Teleseismic  P  waves  recorded  at  EKA  processed  by  beaming  and  adaptive  filtering 
with  and  without  removal  of  site  frequency  responses. 

Figure  46:  Spectra  representative  to  the  signal  beam  loss  (top)  and  the  S/N  gain  obtained  for  two 
events  by  using  the  four  modes  of  processing  described  in  the  text. 
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Figure  47:  Optimum  processing  of  the  nuclear  explosion  Mast,  overlain  by  added  amplified  actual 
background  noise,  at  NORS AR.  The  best  performance,  about  2: 1  S/N  gain,  is  associated  with  site 
equalization  followed  by  adaptive  beaming  (bottom  trace),.  It  must  be  pointed  out,  that  contrary  to 
any  superficial  appearances,  this  is  not  equivalent  to  frequency  filtering  since  this  process  works 
by  spatial  filtering  only  and  does  not  change  the  spectrum  or  waveforms  of  any  signals. 

Figure  48:  Noise  reduction  gain  and  beam  loss  for  western  Norway  earthquakes  as  processed  by 
beaming  and  adaptive  filtering  with  and  without  removing  the  site  transfer  functions. 
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INTRODUCTION 


The  term  "broadband  analysis"  is  often  used  in  seismology  to  designate  time  domain 
studies  of  displacement  seismograms  where  the  effects  of  the  instrument  response  were  made  flat, 
over  a  wide  frequency  band.  The  resulting  seismograms  thus  will  show  the  actual  ground 
displacement  within  a  wide  frequency  band,  but  will  still  contain  the  effects  of  Q  and  source  time 
functions.  Because  of  these  limitations,  matching  the  waveform^  of  such  data  may  not  extract  the 
maximum  information  available,  although  modeling  in  a  wider  frequency  band  is  preferable  to 
working  with  either  the  standard  short  or  long  period  instrument  recordings  alone.  Time  domain 
waveform-matching,  u^mg  RMS  fitting  criteria,  will  inevitably  weight  the  information  in  the 
frequency  domain  proportionately  to  the  actual  power  spectra  of  the  signals  as  they  appear  on  the 
broadband  seismograms.  The  situation  will  be  much  worse  if  we  only  use  the  outputs  of  some 
standard  instruments,  such  as  the  WWSSN  short  period  and  the  LRSM  long  period,  which  have 
notoriously  narrowband  characteristics.  The  problem  is  further  limited  by  the  presence  of 
background  noise.  Assuming  that  we  have  a  knowledge  of  the  variations  of  the  S/N  ratio  as  a 
function  of  frequency  for  the  data  we  wish  to  analyze  the  extraction  of  source  and  path, 
information  from  the  data  can  be  optimized  by  considering  some  basic  formulas  of  information 
theory.  Estimates  of  the  S/N  ratio,  as  a  function  of  frequency,  can  be  obtained  by  simply 
comparing  the  spectra  of  the  seismic  data  in  the  signal  window  to  those  in  the  noise  windows 
preceding  them.  Thus,  we  may  apply  Shannon's  formula  for  the  rate  of  the  information  transfer 
capability  by  a  noisy  channel, 


I  =  BTlog{\+-^) 

*  n 


(1) 


which  basically  tells  us  how  to  weight  the  data  with  respect  to  frequency  in  order  to  extract  the 
most  information.  In  this  equation,  B  is  the  bandwidth,  T  is  the  time  sample  length  and  and 
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are  the  power  levels  of  the  signal  and  noise,  respectively.  If  base  2  logarithm  is  used,  the  formula 
gives  the  amoui  t  of  information  in  bits.  When  the  signal-to-noise  power  ratio  varies  with 
frequency,  we  can  subdivide  the  spectrum  into  bands  with  roughly  constant  S,/N  ratios  and 
estimate  the  information  over  the  total  signal  band  by  summing  the  information  contributions  for  all 
bands.  Weighting  the  data  according  to  the  formula  above,  is  the  ideal  we  should  be  striving  for 
and  it  should  form  the  definition  of  "broadband  processing."  We  should  attempt  to  fit  both  the 
amplitudes  and  phases  cf  the  frequency  components  of  the  signals  with  our  models  throughout  the 
frequency  bands  where  the  S/N  ratio  is  gcod.  It  should  be  clear  that  this  is  not  what  is  being  done 
in  most  seismological  source  inversion  studies.  Waveform  inversion  studies  that  utilize  data  from 
only  one  type  of  narrowband  instrument  and  that  is  uncompensated  for  Q  or  source  effects 
essentially  discard  most  information.  In  effect,  the  information  contained  at  higher  frequencies  is 
often  totally  ignored  and  these  frequencies  often  contain  the  most  information  and  provide  the  most 
resolution  in  space  and  time.  In  the  case  of  misfitted  high  S/N  signals,  we  may  regard  the  power 
in  the  differences  of  the  signals  to  be  fitted  and  those  predicted  from  a  model  as  "misfitting  noise," 
to  be  included  into  Pn  in  Equation  (1). 

The  amount  of  information  that  can  be  extracted  from  conventional  time-domain  source 
inversions  is  determined  by  the  mean  square  effective  bandwidth  fig,  defined  as 

I  f^P,if)df 

^  J.OO 

B,  = -  (2) 

f  p,(f)df 


where  /is  frequency  and  the  integrals  are  centered  around  the  main  peak  (for  unimodal  P(f)).  It 
should  be  clear  that  the  effective  bandwidth  will  depend  on  the  instrument  type,  source  function 
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and  the  Q  along  the  path,  and  thus,  the  results  will  depend  on  these  as  well.  This  is  very 
undesirable,  since  none  of  these  variables  has  anything  to  do  directly  with  the  information  available 
in  the  actual  signals.  If  an  inversion  method  based  on  Equation  (1)  is  applied  and  optimized  over 
various  frequency  bands,  the  results  will  depend  only  on  the  S/N  ratio  and  the  true  information 
content  as  defined  by  Shannon's  formula.  We  shall  return  to  these  points  later. 

An  often  stated  justification  for  the  de-emphasis  of  high  frequency  information  is  that 
waveforms  appear  to  be  quite  unstable  at  frequencies  as  low  as  1  Hz,  and  they  vary  strongly  with 
small  changes  of  location  of  the  sensors  (this  fact,  on  the  other  hand,  did  not  stop  many  workers 
from  fitting,  in  fine  detail,  short  period  waveforms  at  less  than  a  dozen  stations  by  quite  refined 
source-path  models).  It  was  found,  however,  that  these  variations  in  waveforms  can  generally  be 
specified  as  repeatable  path-dependent  inter-sensor  transfer  functions,  which  can  be  estimated  and 
removed  from  the  data  (Filson  and  Frasier,  1972;  Der  et  al,  1985).  Therefore,  it  is  possible  to 
extract  valuable  source  information  from  the  data  at  high  frequencies,  but  at  the  expense  of  much 
more  work  that  has  to  be  done  on  the  estimation  of  site-transfer  functions  and  the  resulting  loss  in 
the  degrees  of  freedom.  High-frequency  seismic  waves  provide  the  most  resolution  with  regard  to 
the  details  of  the  source  processes  of  interest  in  nuclear  monitoring,  such  as  estimating  pP  time 
delays,  strain  release  mechanisms,  etc.,  and  should  not  be  neglected. 

Since  spectral  factorization  is  another  one  of  the  main  themes  of  this  report  and  it  can  be 
applied  to  various  types  of  data,  it  is  necessary  to  include  a  discussion  of  it  prior  to  the  beginning 
of  the  presentation  of  the  results  of  data  analyses.  Spectral  factorability  is  defined  simply  by  the 
following  equation  defining  the  spectra  Fiji co)  of  a  given  short  period  or  high-frequency  seismic 

arrival  from  multiple  events, /■  and  recorded  at  multiple  sensors,;, 


io) )  =  S.  (£y  )Rjio))  +N-J  i(o ) 


(3) 
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where  Si(o})  are  source  spectral  factors,  Rj((o)  are  receiver  (site)  spectral  factors  and  Nij(o))  are 
ambient  noise  terms  (Figure  1).  The  existence  of  site  factors  Rj((o)  is  well-documented  at  seismic 
arrays  and  is  the  main  cause  of  the  variability  of  waveforms  at  the  various  sensors  of  the  arrays. 
The  site  factors  are  the  results  of  the  structural  complexities  near  the  receivers  and  have  been 
shown  to  be  strongly  dependent  on  the  azimuths  and  slownesses  of  arrivals.  Equation  (3)  is, 
therefore,  valid  only  if  the  Rj(03)  are  common  to  the  full  set  of  data  Fij(co),  i.e.,  the  slowness 
vectors  are  the  same  for  all  the  sensors  for  the  data  set.  The  Rj(co)  are,  unfortunately,  not  able  to 
be  modeled  deterministically  from  the  usually  scant  knowledge  of  receiver  structures  and  must  be 
estimated  from  the  data.  It  is  generally  not  possible  to  uniquely  define  both  factors  5/(0))  and 
Rj( co)  from  a  data  set  either,  because  we  can  always  divide  out  some  arbitrary  function,  F( co)  from 
all  5/(0)),  and  multiply  all  Rj(03)  with  it  to  get  new,  equally  valid  source  and  site  functions  (Der  et 
al,  1987).  Nevertheless,  in  some  cases  (teleseismic  P  waves),  some  reasonable  physical  and  semi- 
empirical  assumptions  can  be  made  that  limit  the  possible  choices  and  define  the  two  types  of 
factors  more  uniquely,  thus,  reduci'^g  the  possibilities  of  such  trade-offs.  The  fact  that  Equation 
(3)  holds  in  many  cases  is  hardly  surprising.  One  would  naturally  expect  that  taking  spectral  ratios 
between  waveforms  of  the  same  type  of  pairs  of  events  at  common,  closely-grouped  sites,  would 
tend  to  reflect  the  ratios  of  the  effective  source  functions.  The  main  point  here  is  the  utilization  of 
this  phenomenon  for  the  elimination  of  the  site  effects  and  the  exploitation  of  the  redundancy  in 
such  data  sets  for  reducing  noise.  Clearly,  we  have  redundancy  if  we  have  N>2  events  and  M>2 
sensors,  since  then  MN  >  M+N. 

The  Rj(0)),  besides  being  functions  of  the  azimuth  and  distance  of  the  sources,  will  also  be 
dependent  on  the  type  or  modal  composition  of  the  waves  we  want  to  factor.  This  plays  a  role 
when  we  try  to  analyze  regional  wave  groups,  which  have  common  generic  names  such  as  Lg,  but 
may  consist  of  different  sets  of  modes.  Instead  of  getting  bogged  down  in  the  discussion  of 
various  theoretical  models,  we  shall  attempt  to  clarify  some  of  these  statements  below  with 
concrete  illustrative  cases. 
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The  first  is  the  case  of  P  waves  from  limited  source  regions  recorded  at  individual 
teleseismic  arrays  (Figure  2a).  Since  the  slownesses  and  the  take-off  angles  for  all  arrivals  are 
practically  identical,  the  effective  site-functions  will  remain  the  same  for  the  various  events  at  the 
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FIGURE  1:  Decomposition  and  reconstruction  of  a  matrix  of  seismograms,  using 
site  and  source  spectral  factor  representation. 
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same  sensor.  Moreover,  because  the  array  looks  at  a  very  small  segment  of  the  radiation  patterns, 
the  source  functions  will  be  the  same  for  all  sensors  of  a  given  event.  This  is  the  definition  of 
factorability  in  Equation  (1)  and  thus,  groups  of  all  events,  earthquakes  or  explosions,  will  be 
factorable.  This  is  the  case  which  was  most  extensively  studied  in  the  past  (Der  et  al,  1986;  1987; 
1988).  In  this  application  it  is  assumed  that  the  time  domain  representations  approximately  average 
to  a  delta  function  when  stacked  over  sites.  The  source  functions  thus  derived  will  be  "phaseless 
seismograms"  (Stewart  and  Douglas,  1983)  or  band-limited  representations  of  the  effective  P 
source  function  radiated  in  the  direction  of  the  array. 

Let  us  assume  now  that  instead  of  a  single  array,  we  use  teleseismic  P  wave  recordings  of  a 
number  of  stations  at  roughly  the  same  distance  from  a  limited  source  region,  containing  numerous 
events  (Figure  2b).  Even  if  a  source  has  a  perfectly  symmetrical  pattern  of  radiation,  the 
waveforms  observed  may  vary  considerably  because  the  site  factors  can  be  still  quite  different. 
Nevertheless,  groups  of  events  with  identical  radiation  patterns  will  be  factorable  and  the 
interpretation  of  the  resulting  source  time  functions  remains  valid.  We  have  found  that  the  RSTN 
and  some  subsets  of  UK  array  recordings  of  some  sets  of  Shagan  explosions  were  factorable,  and 
some  groups  of  Degelen  events  were  not  (Der  et  al,  1986;  1987).  On  the  other  hand,  if  we  include 
one  earthquake  into  a  group  of  explosions,  the  group  as  a  whole  will  not  satisfy  the  factorability 
condition  because  of  the  different  symmetry  of  the  radiation  from  earthquakes  making  the  source 
factor  different  for  each  sensor  of  such  events.  The  same  is  the  case  for  the  groups  of  Degelen 
explosions  just  mentioned,  as  seen  at  subsets  of  UK  arrays  which  have  different  radiation 
symmetries,  possibly  because  of  the  complex  topography  at  the  site.  Thus,  factorization  analysis 
in  this  case  can  be  used  to  test  relative  symmetries  in  radiation  of  a  group  of  events  with  obvious 
potential  application  for  discrimination. 

The  third  example  consists  of  recordings  of  a  regional  phases,  such  as  Lg,  from  a  group  of 
events  located  in  a  limited  area  at  a  small  regional  array  (same  situation  as  that  shown  in  Figure  2a). 
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FIGURE  2:  Scenarios  for  which  spectral  factorability  may  be  valid.  a) 
Teleseismic  P  from  a  limited  source  region  recorded  at  a  seismic  array,  b) 
Network  recordings  of  groups  of  events  with  similar  radiation  patterns. 
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The  designation,  "Lg,"  is  based  on  certain  group  velocity  limits,  not  on  the  exact  modal  nature  of 
the  signals.  Consequently,  Lg  wave  groups  can  have  different  modal  compositions  from  sources 
at  similar  locations,  depending  on  the  source  mechanism  and  depths.  Thus,  Lg  wave  groups  with 
similar  modal  compositions  will  interact  with  the  near-site  structures  similarly  and,  thus,  such 
groups  will  be  spectrally  factorable.  In  this  application,  it  would  not  be  prudent  to  assume  that  the 
time  domain  representations  of  site  factors,  averaged  over  sites,  would  still  be  an  impulse.  Thus, 
in  this  case,  we  simply  test  for  the  validity  of  a  factorability  equation  without  attempting  to  interpret 
the  resulting  source  and  site  factors.  When  the  Lg  phases  have  different  modal  compositions,  then 
the  site  factors  that  represent  the  interaction  of  the  different  set  of  modes  and  the  associated  set  of 
site  functions  will  be  different  also.  Thus,  it  can  be  expected  that  groups  of  Lg  {and  other  regional 
phases)  from  events  with  different  modal  compositions,  as  seen  at  the  array,  will  not  satisfy  as  a 
whole  the  factorability  equation  above.  The  obvious  application  of  such  analyses  is  to  group 
closely  located  events  with  respect  to  differences  in  source  mechanisms  or  depth. 

It  is  easy  to  guess  what  the  performance  of  the  process  would  be  in  various  other  cases  of 
source  and  site  configurations  on  the  basis  of  some  general  seismological  principles  and  the 
examples  discussed  above.  Given  that  factorability  is  valid  for  a  given  set  of  data,  the  remaining 
question  is  the  nature  of  the  sets  of  transfer  functions  that  it  implies.  Clearly,  if  Equation  (1)  is 
satisfied  then  it  will  be  possible  to  transform  waveforms  of  the  various  events  at  all  common 
sensors  into  each  other  by  a  filter  constructed  by  taking  the  ratio  of  the  respective  source  functions. 
Similarly,  we  can  transform  the  waveforms  recorded  at  various  sites  into  each  other  for  all 
common  events,  using  filters  that  are  ratios  of  the  respective  site  functions. 

The  factorization  algorithm  used  in  this  report  is  the  same  as  that  used  for  multi-channel 
deconvolution.  A  flow  diagram  of  the  process  is  reproduced  in  Figure  3  for  those  unfamiliar  with 
this  process.  For  further  details  and  statistical  treatments,  the  reader  is  referred  to  our  earlier 
publications  (Shumway  and  Der,  1985;  Der  et  al,  1987).  The  assumptions  regarding  the 
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interpretation  of  source  time  functions  resulting  from  that  procedure  have  no  bearing  with  regard  to 
the  validity  of  factorization  in  the  other  applications  in  Section  I. 

Factorability  implies  that  Equation  (3)  is  satisfied  by  the  spectra  of  a  group  of  events 
(within  limits  allowed  by  the  background  noise).  Ignoring  the  noise  for  the  time  being  (assuming 
high  S/N  ratio),  we  need  some  wav  to  measure  how  well  the  total  energy  in  the  data  can  be 
explained  bv  the  source  and  site  factors  as  a  function  of  frequency.  There  can  be  various  ways  to 
establish  such  a  measure.  Our  choice  is  the  average  ensemble-averaged  coherence  defined  as 

'L'LUr 

C  =  ,  "  '  (4) 

'Z’LUr 

e  s  e  s 

Where  overbar  denotes  spectral  smoothing,  fd  and  /r.are  Fourier  transforms  of  data  and 
reconstructions  and  the  summations  are  over  events  and  sites.  All  quantities  in  Equation  (4)  are 
implicit  functions  of  frequency.  In  order  to  measure  the  trace  similarities  in  small  time 
neighborhoods,  we  apply  heavy  spectral  smoothing.  This  measure  has  the  advantage  that  statistics 
for  testing  significance  are  readily  computed  and  the  effect  of  background  noise  can  also  be 
considered.  Another,  alternative  way  to  measure  the  same  thing  could  be  to  take  the  difference  of 
the  corresponding  data  and  reconstructed  traces  and  compute  the  ratio  of  the  powers  in  these 
residual  traces  and  the  data.  Similar  measures,  with  summations  over  sites  and  events  only,  can  be 
also  used  to  test  the  relative  effectiveness  of  spectral  factorization  for  individual  events  in  a  group 
or  pairs  of  sensors  and  test  alternative  causal  models  of  teleseismic  P  wave  secondary  arrival 
patterns  in  explaining  the  observed  data  (see  Section  2.0). 
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A  data  matrix  of  P  wave  seismogram 
from  a  set  of  events  observed  at 
common  sensors  is  being 
decomposed  into  source  and  site 
factors  as  described  by  the  formula 


Y,j((d)  =  R,(a))  S/oj) 


This  decomposition  is  verified  by 
reconstructing  the  original  data 
from  the  factors  R  and  S. 


FIGURE  3:  Flow  diagram  of  the  spectral  factorization  process  utilized  in  this 
report. 
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1.0  COHERENT  PROCESSING  OF  REGIONAL  SIGNALS  AT  SEISMIC 

ARRAYS 

General  Considerations 

Obtaining  source  information  from  regional  recordings  is  very  important  in  treaty 
monitoring  because  of  the  concern  about  small,  decoupled  nuclear  explosions  possibly  hidden  in 
conventional  firing  patterns  of  quarrying  operations.  Thus  far,  little  progress  has  been  made,  and 
practically  no  work  done,  in  extracting  information  from  regional  waveforms  by  coherent 
processing,  although  spectral  expressions  of  multiple-firing  sequences  have  been  discovered 
(Willis,  1963;  Baumgardt  and  Ziegler,  1988).  The  disadvantage  of  the  spectral  and  cepstral 
methods  is  that  the  phase  information  in  the  signals  is  discarded  and  that  the  modulation  patterns 
depend  on  the  similarity  in  the  waveforms  from  the  various  individual  shots  making  up  the  firing 
sequence.  The  advantage  of  such  methods  is  that  they  are  robust  and  do  not  otherwise  depend  on 
the  details  of  propagation,  unless  multi-pathing  occurs.  In  this  repon  we  shall  discuss  the  general 
problem  of  regional  waveform  analysis  and  show  how  investigating  the  spectral  structure  of  the 
data  may  help  in  recovering  source  information,  in  addition  to  that  obtainable  by  other  analysis 
methods. 

A  general  conclusion  that  may  be  drawn  from  the  numerous  attempts  of  computing 
synthetic  seismograms  of  regional  arrivals,  such  as  Pn,  Pg,  Sn  and  Lg,  and  comparing  them  to 
real  data,  is  that  it  is  quite  impossible  to  compute.  Green’s  functions  with  sufficient  precision  to 
synthetically  reproduce  detailed  features  of  high-frequency  regional  waveforms  use  deterministic 
models.  Waveforms  of  regional  phases  vary  considerably  with  the  locations  of  sensors,  even 
across  small  arrays,  such  as  NORESS.  We  do  not  have  a  sufficiently  detailed  knowledge  of  the 
geological  structures  to  compute  such  waveforms,  and  even  if  we  had  such  knowledge,  we  simply 
do  not  know  how  to  compute  synthetic  seismograms  for  the  complicated  structures  making  up  the 
earth  s  crust  at  the  scale  comparable  to  the  short  wavelengths  involved.  While  it  is  possible  to 
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obtain  synthetic  seismograms  for  long  period  and  some  standard  short  period  seismograms  that 
match  observed  waveforms  in  considerable  detail,  for  regional  seismograms,  the  best  we  can  hope 
for  is  to  match  the  overall  characteristics  of  wave  envelopes  and  relative  phase  amplitudes.  Beyond 
attempting  to  match  such  gross  characteristics  of  regional  arrival,  practically  nothing  has  been  done 
to  utilize  the  information  in  the  details  of  regional  waveforms. 

The  question  arises  whether  it  will  ever  be  possible  to  derive  source  related  information 
from  the  detailed  waveforms  of  regional  seismograms,  if  we  accept  the  fact  that  we  shall  never 
know  the  Green's  functions.  The  answer  to  this  question,  based  on  the  work  to  be  presented 
below,  is  a  qualified  yes,  although  the  possibilities  seem  to  be  limited  due  to  the  extreme  variability 
of  regional  waveforms  with  source  mechanism  and  the  positions  of  the  sources  and  receivers. 

Before  we  present  the  results  of  our  data  analyses,  it  is  necessary  to  describe  some 
expected  spectral  characteristics  of  regional  seismic  arrivals  from  small  explosions  in  the  light  of 
some  elementary  seismological  considerations  and  observational  facts.  We  know,  for  instance, 
that  the  time  durations  of  source  time  functions  of  most  quarry  blasts  are  short,  a  few  tenths  of  a 
second.  In  addition,  spectral  modulations  observed  for  such  events  (Baumgardt  and  Ziegler, 
1987)  could  not  occur  unless  the  subevents  were  very  similar,  so  they  could  coherently  interfere 
with  each  other.  Moreover,  since  the  spectral  modulations  for  all  the  regional  arrivals,  Pn,  Pg,  Sn 
and  Lg,  are  commonly  similar  (Figure  7),  the  spatial  extent  of  such  sources  could  not  be  large. 

Otherwise,  the  modulations  in  the  spectra  would  have  been  different  for  the  various  arrivals 
due  to  the  differences  in  the  delays  in  propagation,  with  the  quite  different  phase  velocities 
associated  with  each  of  these  phases.  Thus,  in  many  cases,  we  may  regard  quarry  blasts  as  point 
sources.  We  observe  that  the  waveforms  of  regional  arrivals  vary  considerably  even  over  a  small 
array.  By  reciprocity,  we  may  expect  that  similarly  small  displacements  of  even  identical  sources 
will  also  cause  differences  in  the  waveforms  if  the  source  region  has  the  same  degree  of  lateral 
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viiriability  in  structures  at  the  observing  array.  In  order  to  utilize  such  information,  the  statistical 
nature  of  these  inter-site  transfer  functions  must  thus  be  studied.  One  would  expect,  similarly  to 
teleseismic  site-transfer  functions,  that  if  the  physical  nature  of  the  arrivals  is  the  same,  such  as  the 
case  for  teleseismic  P,  then  they  would  interact  with  the  geological  structures  near  the  recording 
sites  in  a  similar  fashion.  Thus,  we  would  expect  the  array  spectral  matrices  to  be  factorable  for 
events  in  limited  source  regions  (Filson  and  Frasier,  1972;  Der  et  al,  1987).  The  converse  may  be 
the  case  for  some  regional  arrivals,  such  as  Lg  phases  with  different  modal  composition. 

After  these  preliminary  remarks,  let  us  examine  the  general  spectral  structure  of  seismic 
arrivals  in  an  arbitrary,  heterogeneous  anisotropic  medium  for  various  scenarios.  Let  us  assume 
that  we  have  source  and  receiver  regions  of  limited  spatial  extent  such  that  the  dimensions  of  these 
regions  are  small,  compared  to  the  distance  between  them.  Let  r,-  and  rj  be  the  position  vectors  of 
the  sources  and  the  receivers  in  the  two  regions  (Figure  2).  The  spectra  of  the  observed 
seismograms  from  events  can  be  written  in  the  frequency  domain  by  transcribing  some  known 
equations  (Aki  and  Richards,  1980,  p.  53)  into  the  frequency  domain: 

f „  (® )  =  X  <“  •  r  •  ‘'i )  (5) 

mn 

In  this  equation,  the  are  the  components  of  the  moment  tensor,  each  with  its  own  source  time 
function,  and  the  Gmn  are  derivatives  of  Green’s  functions  (they  will  be  called  Green's  functions 
for  simplicity  in  this  report)  associated  with  them.  The  index,  i,  is  for  the  given  .source  and  j  is  for 
the  sensor.  In  this  expression,  the  index  for  component  in  Aki  and  Richards  (1980),  which  is 
always  vertical  in  this  report,  was  dropped  and  the  summation  over  moment  tensor  components  is 
explicitly  indicated,  rather  than  implied  by  repetition  of  indices,  in  order  to  avoid  confusion  since 
some  other  indices  will  be  repeated  in  the  expressions  in  this  report.  The  di  refer  to  the  dependence 
of  the  Green's  functions  on  .source  depth,  the  receiver  depth  dependence  is  not  indicated  since  the 
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receivers  are  assumed  to  be  on  the  surface.  This  structure  is  equivalent  to  a  multi-channel  filtering 
problem  in  which  we  can  consider  a  seismogram  as  a  result  of  alternatively: 

a)  Filtering  the  source  inputs  for  each  component  of  a  moment  tensor  with 
filters  corresponding  to  the  Green’s  functions. 

b)  Filtering  the  Green's  functions  with  the  six  independent  components  of  the 
moment  tensor. 

It  is  advantageous  sometimes  to  look  at  the  problem  using  the  second  interpretation  since, 
as  we  pointed  out  above,  we  shall  probably  never  know  the  Green's  functions  and  any  information 
we  gain  from  the  seismograms  should  depend  only  on  the  linear  transfer  function  relationships 
among  the  various  seismograms. 

The  relationship  (5)  is  too  general  t  ■  be  of  much  use;  it  cannot  be  factored  in  the  sense  of 
Equation  (3)  unless  something  is  kne  vn  or  a<!<^’!med  about  the  spectral  terms  and  factors  in  this 
expression.  To  see  how  this  equation  can  lead  to  (3)  and  how  Equation  (4)  can  be  utilized  in  some 
applications,  we  must  consider  some  known  facts  and  special  situations.  First  of  all,  we  know  that 
the  moment  tensor  components  generally  have  a  short  finite  impulse  response  (FIR),  relative  to  the 
length  of  the  regional  seismograms,  since  most  seismic  sources  we  consider  are  of  short  duration. 
The  Fourier  transforms  of  short  FIR  functions  will  be  denoted  with  lower  case  in  the  following 
discussion.  Secondly,  we  found  by  experimenting  with  Wiener  filters  of  various  lengths  that  the 
inter-site  transformation  of  waveforms  at  NORESS  could  not  be  efficiently  performed  using  short 
FIR  time  domain  filters.  We  know  that  by  reciprocity  the  displacement  of  sources  of  the  same 
order  as  the  NORESS  array  diameter  will  destroy  waveform  similarity  if  the  source  region  has  the 
same  degree  of  lateral  structural  heterogeneity  as  the  NORESS  site.  We  also  assume  for  the  time 
being,  and  will  show  later  by  data  analyses,  that  the  Green's  functions  factor  into  source  and  site 
factors  for  sources  and  sensors  in  limited  areas.  The  rest  is  a  straightforward  application  of  linear 
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system  theory  to  the  analyses  of  single  and  multi-channel  systems  of  linear  systems  consisting  of 
short  source  function  FIR  filters  (Box  and  Jenkins,  1978)  assumed  to  be  driven  by  "random" 
Green's  functions.  In  the  following,  we  describe  several  scenarios  of  practical  significance.  Due 
to  the  small  inter-sensor  coherence,  we  treat  the  recordings  at  the  various  sensors  as  independent 
realizations,  replications,  of  the  same  random  process. 

We  use  the  term  "coherent  processing"  in  this  report  for  techniques  that  exploit  the 
generalized  factorability  condition  as  defined  above.  This  is  somewhat  different  from  the  usual 
definition  coherent  processing  techniques,  such  as  beaming  and  F-K  analyses,  routinely  applied  to 
regional  signals  at  seismic  arrays  (Mykkeltveit  et  al,  1983;  Ingate  et  al,  1985;  Kwaerna  and 
Mykkeltveit,  1986)  which  basically  depend  on  the  idealized  assumption  of  complete  signal 
similarity.  These  processes  and  the  regional  array  layouts  were  optimized  for  signals  coming  from 
arbitrary  directions  using  the  concept  of  averaged  signal  coherences  vs.  sensor  distances  and  are 
applied  only  at  inter-sensor  distances  and  frequencies  where  the  signal  decorrelation  due  to  site 
effects  can  be  neglected.  However,  if  they  satisfy  equations  of  the  forms  (3)  and  (4),  regional 
signals  across  a  regional  array  are  much  more  "coherent,"  in  a  more  general  sense,  if  we  analyze 
regional  arrivals  from  limited  source  regions.  Let  us  consider  a  few  special  adaptations  of 
Equation  (5)  to  various  practical  situations. 


In  the  case  of  a  pair  of  collocated  events  with  different  source  time  functions,  but  identical 
source  mechanisms  and  depths.  Equation  (5)  can  be  rewritten  as 


Fijico) 


mn  mn  '  J 
mn 


S-{(0) 


(6) 


This  system  can,  thus,  be  modeled  by  two  linear  systems  consisting  of  two  different  source 
time  functions  driven  by  identical  "random"  (Green's  function)  "inputs,"  the  expression  in  the 
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bracket  in  Equation  (6).  '  1*^-  only  difference  between  the  two  sets  of  seismograms  is  that  a 

different  source  time  function  and  the  spectra  are  factorable  in  the  sense  of  Equation  (3),  if  we 
ignore  the  noise  term.  Since  the  source  locations  and  depths  are  assumed  to  be  the  same,  the  Gmn 
will  depend  only  on  rj,  the  receiver  locations,  and  at  each  site  we  shall  have  another  independent 
realization  of  the  same  system.  We  must  note  that,  contrary  to  some  often  expressed  beliefs,  this 
does  not  imply  that  the  waveforms  from  closely-spaced  events  will  have  to  be  similar  or  even  have 
high  average  correlation  coefficients  between  waveforms.  It  is  quite  easy  to  have  essentially  zero 
time  domain  correlation  coefficients  with  a  suitable  selection  of  the  two  source  time  functions. 

A  commonly  observed  phenomenon  for  numerous  quarry  blasts  is  that  the  spectral 
modulation  is  the  same  for  all  regional  arrivals  (Baumgardt  and  Ziegler,  1988),  even  as  seen  at 
several  regional  arrays.  This  situation  corresponds  to  the  situation  described  by  Equation  (6). 
This  implies  a  purely  temporal  modulation,  without  any  spatial  structure  for  the  source  and 
identical  combination  of  Green's  functions.  Without  waveform  similarity  for  the  individual 
subevents,  no  constructive  and  destructive  interference  of  individual  frequency  components  and 
spectral  modulation  could  occur.  Various  scenarios  involving  combined  temporal-spatial  shot 
configurations  of  shots  can  easily  be  simulated  and  such  are  being  actually  applied  for  vibration 
control  in  quarrying  operations  (Anderson  and  Stump,  1989). 

For  this  situation,  the  site  functions  will  be  identical  for  all  the  events,  since  the  modal 
makeup  of  the  wavetypes  for  a  given  regional  phase  for  collocated  events  with  the  same 
mechanism  and  depth  are  expected  to  be  the  same.  The  two  outputs  will  then  be  coherent,  since  it 
will  be  possible  to  transform  one  into  the  other  (at  all  sensors)  with  a  filter  that  is  essentially  si/s2, 
or  its  reciprocal.  Except  for  some  degenerate  cases,  this  filter  will  generally  have  a  short  impulse 
response.  The  existence  of  short  FIR  transfer  function  implies  that  one  will  have  a  high  inter-event 
coherence 
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between  two  events  i  and  j,  in  the  case  of  low  background  noise. 


Central  to  this  idea  is  the  estimation  procedure  for  the  power  spectral  components,  Pij. 
These  must  be  computed  by  some  smoothing  procedure,  such  as  smoothing  the  dot  products  of  the 
signal  Fourier  transforms  denoted  as  vectors  /,  such  that 

Pij  (co)  =  ^fi(co)f*(o))  (8) 

S 

where  the  overbar  denotes  spectral  smoothing,  the  star  denotes  complex  conjugation  and  the 
summing  of  smoothed  products  is  performed  over  sensors.  Smoothing  can  be  accomplished  in 
various  ways,  by  actually  smoothing  the  spectra,  by  averaging  over  shorter  time  windows,  or 
some  combination  of  the  two.  In  our  case,  we  use  the  combination  of  spectral  smoothing  and 
averaging  over  array  sites.  The  equivalent  time-bandwidth  product  (TBWP)  of  the  result  can  be 
computed  by  multiplying  together  the  spectral  averaging  bandwidth  B,  the  total  time  length  of  the 
data  used  for  each  array  site  T  and  the  number  of  array  sites  N.  In  any  coherence  measurement 
some  decision  must  be  made  concerning  the  bandwidth  to  be  used.  High  TBWP  implies  high 
stability  of  the  coherence  results;  with  a  given  amount  of  data  this  can  be  traded-off  against  the 
resolution  in  spectral  detail  (Bendat  and  Piersol,  1966;  Shumway,  1988).  In  the  case  discussed 
here,  we  can  assume  that  the  waveform  differences  are  associated  with  different  effective  source 
functions,  i.e.,  finng  sequences.  As  we  mentioned  before,  these  generally  have  very  short  time 
durations.  The  amount  of  smoothing  that  can  be  applied  in  coherence  calculations  depends  on  the 
time  length  of  the  impulse  response  of  the  transfer  function  between  the  two  time-series.  As  a  rule 
of  thumb,  one  can  smooth  over  a  bandwidth  of  7/7/  where  7/  is  the  duration  of  such  an  impulse 


response.  This  implies  that  heavy  spectral  smoothing  may  be  applied  in  the  calculation  of  C  above 
if  the  differences  in  waveforms  are  mostly  due  to  different  firing  sequences.  For  further 
discussions  of  the  particulars  of  coherence  calculations,  we  refer  the  reader  to  Shumway  (1988, 
pp.  70-73). 


The  second  case  we  consider  is  that  of  two  or  more  events  with  identical  source 
mechanisms  and  depths,  but  with  larger  differences  in  source  location.  In  this  case,  the  spectra  of 
a  regional  arrival  from  these  events  may  be  still  factorable,  as  described  by  the  equation 


Fijio)) 
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S,  (CO ,  r . )  s-  (CO ) 


(9) 


but  any  pair  of  events  is  no  longer  coherent  in  the  sense  of  Equations  (7)  and  (8).  The  inter-event 
transfer  function  is  no  longer  characterized  by  a  short  FIR  filter  because  of  the  factor  Sjf  co)  now 
implies  the  longer  FIR  filter  due  to  source  dislocation  (a  reciprocal  equivalent  of  array  site  effects). 

In  the  case  of  collocated  events  with  different  source  mechanisms.  Equation  (5)  can  be 
rewritten  in  yet  another  form 


(0))  =  ^  (Q) )  (CO ,  r,  )  (10) 

mn 

In  this  case,  the  waveforms  of  individual  events  will  not  be  coherent  in  the  sense  of  Equation  (6) 
and  (10)  is  not  factorable  either.  Moreover,  the  spectra  of  the  superpositions  of  individual  shots 
with  different  mechanisms  will  not  be  modulated  in  any  easily  predictable  fashion.  A  possible 
scenario  of  such  a  situation  in  quarry  blasting  is  a  mix  of  shots,  some  of  which  may  be  applied  to  a 
vertical  rock  face  (shear  source)  and  some  buried  under  a  horizontal  surface,  both  kinds  of  shots 
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being  close  to  the  surface  at  roughly  the  same  depth  but  involving  the  various  Gij  (co)  with  different 
weights. 

Note  that  the  components  of  the  moment  tensors  in  this  equation  are  written  with  a  lower 
case  to  designate  the  fact  that  they  are  short  FIR  filters.  This  is  a  multi-channel  filtering  scenario 
and,  as  long  as  the  six  independent  Green’s  functions  are  common  to  a  set  of  any  seven  events, 
they  will  be  coherent  in  the  multiple  coherence  sense,  i.e.,  we  shall  be  able  to  reconstruct  the 
waveform  of  the  i-th  time  series  from  the  rest  if  all  have  distinctly  different  source  mechanisms 
(Bendat  and  Piersol,  1966). 

We  believe  that  this  is  a  worst  case  scenario,  since  not  all  six  independent  G  will  contribute 
significantly  to  the  waveforms  and  in  many  cases  only  a  few  may  be  important.  In  such  cases, 
multiple  coherences  with  lesser  number  of  events  may  attain  relatively  high  values  that  are 
significantly  different  from  zero.  Thus,  no  matter  how  complex  the  source  mechanisms  may  be,  it 
will  be  possible  to  "calibrate"  parts  of  a  quarry  with  a  small  number  of  events  and  the  high  multiple 
coherences  will  identify  all  subsequent  events  which  are  roughly  at  the  same  depth  and  location 
within  that  quarry.  In  this  study  we  shall  not  present  any  multiple  coherence  results. 

The  last  case  is  when  the  events  are  too  far  from  each  other  so  that  the  site  effects  at  the 
same  sensors  are  all  different  and  we  are  back  to  the  general  Equation  (5)  again.  This  case  is 
analogous  of  testing  outputs  of  several  different  systems  with  independent  random  inputs  for 
coherence.  These  should  test  as  not  being  significantly  different  from  zero  and  nothing  more  can 
be  said  about  these  events.  This  case  also  covers  all  events  with  different  depths.  If  all  Gmn  are 
different,  then  none  of  the  single  and  multiple  coherences  between  that  event  and  the  events  at  a 
"calibrated"  quarry  will  be  high. 
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The  last  case  corresponds  to  the  practically  important  situation  from  the  monitoring  point  of 
view  when  an  attempt  is  made  to  hide  a  decoupled,  larger  explosion  among  normal  quarry  blasts, 
presumably  at  significantly  different  source  depths,  in  a  quarry.  In  this  case,  the  lack  of  single  or 
multiple  coherence  and  factorability  could  be  used  as  a  criterion  for  detecting  such  suspicious 
events.  Given  the  fact  that  decoupled,  non-chemical  explosions  have  to  be  buried  deeper  than 
conventional  explosives  in  standard  quarrying  procedures,  and  that  the  effective  mechanisms  for 
the  two  types  of  events  will  be  no  doubt  different,  it  does  not  seem  to  be  very  easy  to  mix  the  two 
types  of  explosions  in  order  to  hide  testing. 

The  advantage  of  applying  simple  linear  system  analysis  methods  is  that  one  does  not  need 
to  know  the  Green's  functions  or  the  moment  tensor  time  functions,  which  are  next  to  impossible 
to  compute  with  sufficient  realism  for  regional  arrivals.  We  simply  test  the  multi-sensor  data, 
treated  as  independent  realizations  of  the  same  processes,  for  the  existence  of  linear  transfer 
function  relationships  of  various  types  described  above.  The  types  of  relationships  found  can  be 
related  to  various  relative  source  mechanism  and  location  scenarios.  This  approach  has  the 
advantage  that  it  can  employ  some  simple,  standard  time-series  analysis  techniques  and,  unlike 
most  methods  applied  to  regional  arrivals,  exploit  detailed  regional  seismic  waveform  information. 
Theoretical  developments,  conceptualizations  and  observational  results  of  a  similar  nature,  based 
on  linear  system  theory,  were  also  presented  in  a  forthcoming  report  by  Harris  (1989). 

A  summary  of  the  diagnostics  for  grouping  of  events  with  respect  to  location 

mechanism  and/or  source  depth 

A  short  summary  of  the  expected  features  of  various  scenarios  discussed  above  is  as 
follows: 

a)  Events  with  very  similar  or  identical  waveforms  for  all  regional  arrivals  - 
collocated  events  with  closely  identical  source  time  functions. 
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b)  Factorable,  inter-event  transfer  functions  can  be  characterized  by  relatively 
short  FIR  wavelet  shaping  filters  (Robinson,  1967)  and  these  are  valid  for 
several  regional  phases,  inter-event  coherences  are  high  -  closely-spaced 
events  with  similar  source  mechanisms  and  source  depths,  but  different 
source  time  functions. 

c)  Factorable,  but  short,  FIR  filters  cannot  describe  the  inter-event  transfer 
functions  (low  inter-event  coherences)  -  events  with  similar  mechanisms 
and  depths,  still  grouped  tightly,  but  not  very  closely  spaced. 

d)  Not  factorable,  but  high  multiple  coherences  within  a  larger  group  of  events 
-  closely  spaced  events  with  different  mechanisms. 

e)  Not  factorable,  low,  multiple  coherences  -  any  combination  of  different 
depths  and/or  large  spatial  separations. 


Analyses  of  NORESS  Data  from  Quarry  Blasts  in  Scandinavia 


Coherence  and  cross-filtering  analyses 

In  the  following,  we  shall  demonstrate  by  analyses  of  actual  regional  phases  that  many  of 
the  phenomena  and  properties  predicted  above  are  actually  exhibited  by  regional  data.  Properties, 
such  as  spectral  factorability  for  regional  arrivals,  cannot  be  taken  for  granted  just  on  the  basis  of 
theoretical  speculations  such  as  those  presented  above.  We  were  pleasantly  surprised  to  find  that 
many  of  our  expectations  were  confirmed  despite  the  extremely  complex  nature  of  regional  phases 
and  the  crustal  structures  through  which  they  propagate. 


Our  data  came  from  NORESS  recordings  of  events,  listed  in  Tables  1  and  II,  identified  as 
quarry  blasts  at  the  Titania  mine  and  various  Estonian  mines  named  El,  E4  and  E9.  We  have  also 
looked  at  one  event  at  the  Blasjo  dam  site.  We  do  not  know  about  the  location  accuracy  of  these 
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events  or  the  reasons,  other  than  location,  for  their  association  with  any  particular  mine.  The 
general  locations  of  the  events  analyzed  are  shown  in  Figure  4. 

In  order  to  apply  any  of  the  ideas  above,  we  must  first  verify  the  applicability  of  the 
factorability  condition  for  any  data  set.  We  have  applied  the  multi-channel  factorization  technique 
described  by  Shumway  and  Der  (1985)  and  Der  et  al  (1987;  1990)  to  two  data  sets,  the  Pn  arrivals 
from  this  group  of  Titania  mine  blasts  and  the  Lg  arrivals  from  a  set  of  mine  blasts  at  various 
Estonian  mines.  The  decomposition  of  data  into  source  and  site  factors  is,  of  course,  not  unique 
(see  Der  et  al,  1987),  but  this  is  not  a  major  concern  here.  We  merely  want  to  test  for  the  existence 
of  factorability,  one  of  our  diagnostics  for  source  scenarios,  by  reconstructing  the  complete  set  of 
original  waveforms  from  the  much  smaller  set  of  estimated  source  and  site  spectral  factors  and 
comparing  the  result  to  the  original  data.  High  similarity  verifies  that  the  data  has  the  spectral 
structure  described  by  Equation  (3).  Figures  5  and  6  show  some  representative  time  domain 
examples  of  these  reconstructions  for  Pn  from  quarry  blasts  at  the  Titania  mine  and  Lg  phases  from 
the  Estonian  E4  mine,  respectively.  All  of  these  have  similar  qualities  within  each  set  (showing 
whole  data  sets  would  simply  not  be  practical  because  of  the  large  number  of  traces  involved). 
Typically,  the  correlation  coefficients  between  the  pairs  of  original  and  reconstructed  traces  were 
near  or  above  0.8,  which  is  not  as  good  as  in  the  teleseismic  cases  for  P  waves,  but  still 
respectable,  considering  the  high  signal  frequencies,  up  to  15  Hz,  we  utilized  here  and  the  high 
degree  of  crustal  heterogeneity  at  the  wavelengths  involved  and  the  fact  that  the  recordings  were 
not  noise-free.  Note  that  the  waveforms  are  quite  different  at  the  various  sensors  for  the  same 
event  and  for  the  various  events  at  the  same  sensors.  Joint  factorization  of  the  Titania  events  and 
the  Blasjo  event  did  not  result  in  acceptable  reconstructions,  showing  that  an  azimuth  difference  of 
10  is  too  large  for  using  the  same  site  factors.  In  all  the  work  we  present  in  this  report,  we  have 
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TABLE  I:  List  of  “Titania”  Events 


Event 

Date 

Oriain  Time 

Latitude 

Longitude 

Magnitude 

Titania  1 

14  Feb  86 

14:13:24.9 

58.3 

6.4 

2.7 

Titania  2 

14  Feb  86* 

17:54:10.6 

58.3 

6.4 

2.3 

Titania  3 

07  Jan  86 

14:14:28.9 

58.3 

6.4 

2.2 

Titania  4 

17  Jan  86# 

14:11:01.5 

58.3 

6.4 

2.7 

*  Probably  removed  from  or  not  at  the  same  part  of  quarry. 

#  Reference  event  (unmodulated). 


TABLE  II:  List  of  “Estonia”  Events 


Event 

Date 

Origin  Time 

Latitude 

Longitude 

Magnitude 

E4-1 

19  Jan  88 

11:55:35.0 

59.3 

27.2 

2.5 

E4-2 

20  Jan  88 

12:23:06.0 

59.3 

27.2 

2.5 

E4-3 

21  Jan  88 

11:25:30.0 

59.3 

27.2 

2.4 

E4-4 

26  Jan  88 

09:55:37.0 

59.3 

27.2 

2.4 

EM 

19  Jan  88 

11:50:50.0 

59.3 

24.4 

2.5 

El-2 

21  Jan  88 

12:53:23.0 

59.3 

24.4 

2.2 
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Location  of  the  quarries 


FIGURE  4:  Map  of  Scandinavia  with  the  location  of  NORESS  and  Titania,  Blasjo 
and  the  Estonian  mines  studied. 
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FIGURE  6:  a)  Selected  reconstructed  and  data  traces  of  Lg  arrivals  from  the  E4-1 
event,  b)  The  same  for  the  E4-4  event.  The  designation  C2  refers  to  the  second 
sensor  in  the  C  ring  of  NORESS,  etc. 
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used  12  sensors  of  NORESS  chosen  such  that  they  were  roughly  evenly  disri.buted  over  the  area 
of  the  itrray. 

In  Figure  18,  we  show  some  factoring  results  in  terms  of  coherences  for  a  group  of 
FENNOLORA  shots  at  NORSAR  subarray  3C.  The  efficiency  of  factoring  at  this  subaTay  with  5 
km  sensor  spacing  indicates  that  regional  arrays  larger  than  NORESS  may  be  effective  when  the 
more  general  concept  of  signal  coherence,  implied  by  Equation  (3),  is  applied.  For  instance,  site- 
equalization  filtering  could  be  applied  at  such  arrays  prior  to  beaming  or  optimum  filtering,  thus 
leading  to  better  estimation  of  source  characteristics. 

Starting  with  the  Pn  data  from  the  Titania  mine,  let  us  assume  initially  that  all  these  events 
had  the  same  source  mechanism  and/or  location  and  test  the  idea  that  the  differences  in  waveforms 
are  totally  attributable  to  differences  in  scarce  functions  co)  implied  in  the  firing  sequences 
(Equation  6).  Inspecting  the  spectra  of  mine  blasts  from  the  Titania  mine  in  Norway  at  NORESS, 
it  can  be  observed  that  regional  phases  from  some  of  these  events  are  strongly  modulated,  while 
spectra  of  some  others  are  quite  smooth  (Figure  7).  These  modulation  features  are  similar  for  all 
phases,  Pn,  Pg,  Sn  and  Lg.  This  points  to  a  common,  temporal  modulation  for  each  event,  rather 
than  to  a  pattern  imposed  by  the  spatial  configuration  of  charges,  since  in  the  latter  case  we  should 
see  different  modulation  patterns  for  the  various  regional  phases  which  possess  notably  different 
phase  velocities.  Although  the  lack  of  modulation  does  not  rule  out  some  very  complex  spatio- 
temporal  modulation  patterns  of  the  source  (although  these  must  be  quite  unlikely),  such  that  the 
spectral  manifestations  of  it  are  suppressed,  it  is  reasonable  to  assume  that  the  spectrally 
unmodulated  quarry  blasts  have  simple  impulsive  source  time  functions  or  very  small  delays  while 
the  modulated  ones  have  been  ripple-fired. 

Let  us  compute  now  the  inter-event  coherences  according  to  equations  (5)  and  (6)  above, 
summing  the  smoothed  auto-  and  cross-spectral  estimates  over  sites  only,  for  Pn  phases  from 
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several  pairs  of  these  events  including  the  relatively  unmodulated  one.  The  results  show  (Figure  8) 
that  the  pairs  Titania4-1  and  Titania4-3  have  relatively  high  average  coherences,  indicating  nearly 
identical  source  mechanism  and/or  location,  while  the  pair  Titania4-2  has  lower  coherence, 
indicating  that  event  TitaniaS  has  probably  a  different  mechanism  or  location.  Although  the  S/N 
ratio  was  the  also  the  lowest  for  Event  2,  this  coherence  still  seems  to  be  still  too  low  for  a  coherent 
signal.  The  TBWP  for  these  coherences  was  at  least  120. 

Let  us  hypothesize  now  that  the  unmodulated  event  Titania4  is  a  simple  one,  a  single  shot. 
Consequently,  the  time  domain  impulse  response  of  the  moving  average  (MA)  filter  that  can 
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FIGURE  7:  Spectra  for  various  regional  arrivals  from  the  Titania  mine, 
that  the  spectral  modulation  is  similar,  regardless  of  the  type  of  arrival. 
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transform  the  waveforms  of  this  event  to  those  of  the  Titanial  must  be  a  band-limited 
representation  of  the  firing  sequence.  Therefore,  we  designed  a  Wiener  filter  to  transform  Pn 
waveforms  of  event  Titania4  into  those  of  event  Titanial  at  all  sensors  of  NORESS.  Prior  to 
applying  this  technique,  we  wanted  to  ensure  that  the  traces  processed  are  as  free  as  possible  from 
background  noise  and  other  distortions  of  the  spectrum  by  applying  a  common,  minimum  phase 
band-pass  filter  to  all  traces  of  both  events  that  emphasized  the  energy  and  flattened  the  spectra  in 
the  3  -  15  Hz  frequency  band  where  the  signal-to- noise  ratio  was  higher. 

Computation  of  time-domain  MA  type  filters  can  be  easily  accomplished  by  several  well- 
known  algorithms,  using  computer  codes  widely  available  (Marple,  1987).  Time  domain  design  is 
convenient  because  we  want  to  limit  the  length  of  the  filter.  The  Pn  arrivals  were  lined  up  for  both 
events,  windowed  and  tapered,  and  the  auto-  and  cross-correlation  functions  needed  by  the  filter 
design  algorithms  were  computed  by  ensemble-averaging  these  correlation  functions  over  sites. 
The  ensemble-averaging  process  of  correlation  functions  over  common  sites  suppresses  the  effects 
of  the  site  functions  and  much  reduces  those  of  the  ambient  noise.  The  length  of  the  portion  of  the 
cross-correlation  function  with  highest  amplitude  gave  a  good  indication  of  the  maximum  length  of 
the  ripple-firing  sequence,  since  the  time  length  of  the  large  amplitudes  in  this  function  cannot 
exceed  this  length.  For  the  cases  discussed,  this  upper  limit  is  roughly  of  the  order  of  .5  seconds. 
The  compound  correlation  functions  constitute  the  inputs  to  the  Levinson  algorithm  used  for  the 
computation  of  the  cross-equalizing  Wiener  filter,  consisting  of  3 1  weights. 

Subsequently,  the  filter  was  applied  to  the  simple  event  to  derive  waveforms  of  the 
complex  events  at  all  sites.  A  set  of  traces  showing  the  representative  improvement  between  the 
two  events  at  common  sensors  are  shown  in  Figure  9.  [Note  that  the  original  waveforms  are  much 
different  for  the  two  events  and  waveforms  became  more  similar.]  Nevertheless,  the  increase  in 
the  time  domain  correlation  coefficients  was  only  modest,  from  0.59  to  0.75  for  the  first  pair,  and 
0.49  to  0.57  for  the  second.  This  indicates  that  either  we  need  a  longer  filter  or  that  the  event  pairs 
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FIGURE  9:  a)  Wiener  filtering  results  for  transforming  Event  Titania4  into  Event 
Titanial. 
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FIGURE  9:  b)  Wiener  filtering  results  for  transforming  Event  Titania4  into  Event 
Titania3. 
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are  too  far  from  each  other  and,  thus,  we  cannot  attribute  the  waveform  difference  entirely  to  the 
source  time  functions.  The  same  method  did  not  work  at  all  for  the  pair  Titania4-2. 

We  have  also  applied  the  cross-equalization  for  the  whole  wavetmins,  consisting  of  all  the 
regional  arrivals  of  events  Titania4-1.  The  processing  resulted  in  a  spectacular  increase  of  the 
correlation  coefficients  between  the  two  events  from  .08  to  .65  over  the  whole  set  of  sensors, 
again  confirming  that  the  waveforms  of  these  events  are  closely  related. 

Having  established  that  such  techniques  may  be  suitable  to  identify  groups  of  events  with 
nearly  identical  mechanisms  and/or  locations,  we  examined  the  nature  of  the  relative  inter-site 
transfer  functions  for  the  same  event.  We  know  from  spectral  factorization  that  such  transfer 
functions  exist,  and  they  are  consistent  event  to  event.  Otherwise,  we  would  not  have  been  able  to 
jointly  deconvolve  and  reconstruct  all  these  events.  It  would  be  desirable  if  these  transfer  functions 
would  be  simple,  similarly  to  the  inter-event  transfer  functions  for  the  coherent  event  pair,  because 
then  relative  displacement  over  distances  of  the  order  of  the  NORESS  diameter  would  not  destroy 
the  inter-event  coherence  (applying  the  principles  of  reciprocity).  In  that  case,  the  lack  of 
cohtic.ice  among  events  in  a  limited  source  region  observed  at  an  array  would  be  mostly  due  to 
differences  in  source  mechanisms. 

Unfortunately,  this  is  not  the  case.  We  have  tried  to  design  short,  31  point,  time  domain 
filters  for  equalizing  Pn  waveforms  at  sensors  D1  and  D5  for  the  suite  of  four  events,  but  they 
essentially  failed  to  produce  even  remotely  similar  waveforms.  Thus,  the  inter-site  transfer 
functions  are  generally  complex,  even  over  the  small  NORESS  diameter,  with  very  long  time 
domain  impulse  respon.ses.  Applying  the  reciprocity  argument,  this  implies  that  a  comparably 
small  displacement  in  the  vocation  of  a  source  would  destroy  the  waveform  similarity,  assuming 
that  the  source  region  and  the  .NORESS  site  are  similar  in  the  degree  of  structural  heterogeneity. 
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Site  equalization  processing 


A  major  problem  in  locating  regional  events  is  that  the  first  arrivals  are  often  small  and 
buried  in  noise  and  that  most  phases  are  emergent,  without  clearly  defined  arrival  times.  Azimuths 
estimated  from  F-K  analyses  are  often  not  very  accurate,  since  the  site  effects  tend  to  destroy  the 
plane  wave  character  of  the  signals,  thus  broadening  the  main  lobes  of  array  response  patterns. 
There  are  conflicting  demands  of  keeping  the  array  apertures  small  for  ensuring  signal  similarity 
and  increasing  the  directional  resolution  of  the  arrays.  The  fact  that  the  site  transfer  functions  are 
consistent,  albeit  not  simple,  still  leaves  open  the  possibility  of  utilizing  this  internal  consistency 
for  refining  the  estimation  of  relative  azimuths  among  closely  spaced  events.  If  the  forms  of  the 
inter-site  transfer  functions  do  not  change  much  with  small  azimuth  changes,  except  for  some  small 
time  shifts,  then  such  a  property  could  be  exploited.  To  lest  this  idea  we  have  band-pass  filtered 
the  four  Titania  events  to  emphasize  the  band  with  the  maximum  signal-to-noise  ratio,  between  3- 
10  Hz.  Subsequently,  we  have  taken  the  frequency  domain  site  factors  derived  from  the  multi¬ 
channel  deconvolution  calculations  and  multiplied  the  Fourier  transform  of  trace  D1  with  the  site 
factor  of  D5  and  vice  versa.  This  is  basically  an  intercorrelation  technique  applied  to  the  broadband 
site  (rather  than  event)  equalization.  The  resulting  spectra  were  then  correlated  in  the  frequency 
domain  and  transformed  back  into  the  time  domain 
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to  view  the  cross-correlation  function.  In  this  expression,  the  d's  are  the  Fourier  transforms  of  the 
traces  and  the  R's  are  the  site  transfer  functions.  The  results  for  three  of  the  events  are  shown  in 
Figures  7  and  8.  Event  2  was  too  noisy.  In  inspecting  these  figures,  we  must  point  out  that  the 
results  of  cross-equalization  do  not  have  to  look  like  normal  seismograms,  with  clearly  defined 
phases.  In  multiplying  with  the  complex  site  factors  we  essentially  perform  a  circular  convolution 
with  considerable  wraparound  effect.  It  is  interesting  to  note  that  the  correlation  peak  for  event 
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Titania4  is  not  lined  up  with  those  of  Titania3  and  Titanial  (Figure  10),  indicating  the  possibility  of 
a  small  time  shift.  Figure  1 1  demonstrates  that  the  site-equalized  trace  waveforms  indeed  became 
quite  similar.  The  shift  in  the  peak  does  not  seem  to  be  associated  with  noise  or  mismatch,  since 
all  cross-correlation  functions  for  this  set  of  events  have  a  well-defined,  single  maximum  with  litde 
ambiguity  in  value  of  the  time  shift. 

The  same  process  was  repeated  with  the  recordings  of  Lg  phases  from  quarry  blasts  in 
Estonia  as  recorded  at  NORESS.  The  Lg  phases  of  the  Estonia  mine  blasts  are  somewhat  noisy 
and  the  S/N  ratio  is  only  good  below  6  Hz;  the  rest  of  the  phases  are  buried  deeper  in  noise  than 
those  for  the  Titania  mines.  Applying  the  process  to  four  E4  mine  events,  we  obtained  good  trace 
equalizations,  with  all  correlation  peaks  lined  up  at  zero  lag  (Figure  12).  Using  four  events  from 
E4  and  two  from  El  mines  in  a  joint  .iite-equalization  process,  we  have  found  that  the  efficiency  of 
the  inter-site  equalization  drastically  deteriorated,  as  indicated  by  the  decrease  in  the  correlation 
coefficients  between  equalized  traces  (Figure  13,  left  side).  The  correlation  peaks  are  still  generally 
at  the  same  times,  but  one  of  the  peaks  (for  event  El-2)  is  shifted  in  time.  The  decrease  in  site- 
equalization  efficiency  may  be  attributed  to  the  fact  that  mines  E4  and  El  are  too  far  from  each 
other  to  have  identical  site  functions  at  NORESS.  The  waveform  similarity  is  not  even  visible  for 
this  set  (Figure  1 3,  right  side). 

Given  the  fact  that  after  applying  some  cross-equalization  treatment  to  the  outputs  of  a  pair 
of  sensors  located  at  the  extremities  of  a  small  array  for  events  in  some  limited  source  areas  (not  the 
two  Estonia  mines  together),  the  waveforms  become  quite  similar,  it  may  not  be  justifiable  to  limit 
the  sizes  of  regional  arrays  to  a  few  kilometers.  It  appears  that  regional  signals  from  limited  source 
regions  are  actually  "coherent,"  in  the  sense  of  a  broader  definition  of  this  term,  over  much  larger 
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FIGURE  10:  Cross-correlations  of  the  site-equalized  traces  of  Pn 
sensors  D1  and  D5  for  selected  Titania  events. 
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FIGURE  11:  Site-equalized  Pn  traces  for  selected  Titania  events  between  sensors 
D1  and  D5. 
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FIGURE  12:  Cross-correlation  of  the  site-equalized  traces  and  the  equalized 
traces  of  Lg  between  sensors  D1  and  D5  for  the  E9  Event.  Note  the  excellent 
similarity  between  the  cross-equalized  traces. 
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FIGURE  13:  Results  of  site-equalization  when  spectral  factors  from  joint 
factorization  of  El  and  E4  Events  were  utilized.  Although  the  cross-correlations 
still  have  peaks  at  the  proper  alignment  with  reduced  correlation  coefficients,  (a) 
the  efficiency  of  site-equalization  is  greatly  degraded  as  evidenced  by  the  poor 
similarity  of  the  site-equalized  traces  (b). 
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distances  between  sensors.  Clearly,  the  maximum  array  diameters  dictated  by  the  average 
waveform  correlation  properties  (Mykkeltveit  et  al,  1983;  Der  et  al,  1988)  do  not  applv  here.  It 
could  be  fruitful  to  test  and  develop  similar  techniques  using  data  from  small  arrays  larger  than 
NORESS. 

Factorizability  studies  and  coherence  grouping  of  Estonian  mine  blasts 

Thus  far,  we  have  only  presented  various  time  domain  results  using  inter-event  and  inter¬ 
site  filtering  and  factorization,  demonstrating  that  there  exists  a  remarkable  degree  of  coherency  in 
the  detailed  waveforms  of  regional  Pg  and  Lg  arrivals.  In  the  following,  we  test  the  effectiveness 
of  spectral  factorization  for  Lg  phases  from  Estonian  mine  blasts  by  measuring  the  efficiency  of 
reconstructions  of  data  from  the  estimated  site  and  event  factors  for  various  groupings  of  these 
events.  After  factoring  and  reconstructing  the  waveforms  of  various  groups  of  events,  we 
compute  the  coherence  measure  (4)  between  the  data  and  reconstructions  for  each  event,  using  the 
data  for  all  sensors  in  the  2-6  Hz  band  where  the  S/N  ratio  was  the  best  summing  only  over  sites. 
The  S/N  ratios  for  the  chosen  events  were  such  that,  in  this  frequency  band,  (4)  should  have 
resulted  in  coherences  of  at  least  0.9,  had  the  factoring  on  the  signal  part  been  effective  for  all 
events.  The  statistic  measures  the  similarity  of  the  reconstructions  vs.  the  original  data.  We  are 
not  interested  in  either  the  source  or  site  factors  themselves,  which  we  have  no  way  to  interpret, 
only  in  the  efficiency  of  the  factorability.  It  is  expected  that  events  with  identical  site  factors 
(similar  Lg  modal  compositions  -  hence  similar  mechanisms  and  depths)  will  factor  efficiently. 
The  converse  is  also  true,  events  with  markedly  different  mechanisms,  depth  and  location  will  not 
factor  well  when  grouped  with  those  that  are  close  in  these  parameters. 

We  found  examples  that  confirm  types  of  behavior  that  we  expect  on  physical-theoretical 
grounds.  In  Figure  14,  we  show  a  .set  of  coherences  computed  with  equation  (4)  summing  over 
sites  only  for  any  given  event,  for  four  events,  three  at  the  E9  mine  and  one  at  the  E4  mine.  The 
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FIGURE  14:  Site-averaged  coherence  results  from  the  joint  factorization  of  three 
E9  events  and  one  E4  event. 
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E4  mine  is  only  about  0.2  degrees  off  in  azimuth  at  NORESS,  but  is  1 10  km  distant  from  the  E9 
mine.  The  coherence  for  the  E4  mine,  when  factored  with  the  three  E9  events,  is  visibly  less  than 
for  the  E9  group.  Apparently,  the  site  effects  for  this  event  are  different,  possibly  because  of  the 
differences  in  the  nature  of  Lg  wavetrain.  This  example  shows  that,  in  spite  of  the  very  small 
difference  in  azimuth  that  could  not  possibly  have  been  resolved  from  the  Lg  by  the  NORESS 
array  (Harris,  1990),  we  can  show  that  something  is  different  about  this  event  utilizing  its  Lg 
phase  alone.  This  case  has  limited  practical  significance,  since  we  could  locate  this  blast  at  a 
different  mine  from  other  sources. 

A  more  interesting  observation  is  that  not  all  Lg  phases  from  quarry  blasts  that  were  listed 
to  have  been  fired  at  the  E9  mine  factor  efficiently  when  grouped  together.  We  have  found  two 
factorable  groups.  Figure  15  shows  the  factorization  efficiency  results  for  three  E9  events  from 
one  of  the  groups  (on  the  right)  factored  with  one  from  the  other  group  (left).  The  figure  shows 
that  there  is  something  different  about  this  event  on  the  left  and  that  the  coherence,  instead  of  being 
near  0.9  as  expected  based  on  the  prior  estimate  of  the  S/N  ratio,  is  closer  to  0.6.  In  order  to  find 
what  is  different  between  the  two  groups  we  went  back,  after  the  fact  of  their  discovery,  to 
examine  the  bandpass-filtered  trace  envelopes  for  all  the  E9  events.  We  have  found  that  with 
regard  to  the  envelope  shapes  of  Lg,  they  also  fall  into  two  subtly  different  groups  (Figures  16  and 
17).  For  Type  1  (corresponding  to  coherence  diagrams  on  the  right  in  Figures  16  and  17),  the  Lg 
reaches  a  maximum  suddenly  and  decreases  gradually,  while  for  Type  2,  it  has  a  gradual  beginning 
and  reaches  the  maximum  later.  The  timing  of  the  discernable  phases  does  not  seem  to  be 
different,  thus  agreeing  with  the  idea  that  they  were  at  the  same  mine.  The  first  impression  was 
that  it  may  be  a  difference  in  the  signal-to-noise  ratios  that  could  cause  the  differences  in  the 
coherences.  Indeed,  the  Type  I  as  a  group  is  less  noisy.  But  if  we  compare  the  S/N  in  Lg  for 
events  in  Type  II  and  E95,  they  are  not  much  different.  This  points  to  the  alternative  interpretation 
that  the  modal  compositions  of  the  Lg  phases  in  the  two  groups  were  different  (and  thus,  their 
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mechanisms).  The  other  possibility  that  the  events  were  mislocated  cannot  be  ruled  out  either, 
since  this  occurs  occasionally. 

The  question  of  optimum  array  size  for  regional  events  was  touched  on  previously.  An 
interesting  observation  concerning  the  factorization  of  Pn  phases  is  shown  in  Figure  18.  These  are 
coherence  results  (site-  and  event-averaged)  for  a  group  of  closely  spaced  FENNOLORA  profiling 
shots  at  the  3C  subarray  of  NORSAR.  Given  the  large  spacing  of  these  sensors  (~5  km),  this 
shows  that  the  waveforms  of  these  events  are  "coherent"  in  a  wider  sense. 

These  results  show  that  simple  joint  analyses  of  waveforms  of  regional  arrivals  from 
groups  of  events  can  potentially  be  quite  useful  for  ascertaining  what  activities  take  place  within 
small  source  regions  and  that  evasion  attempts  may  be  detectable.  What  is  needed  in  the  future  is 
to  obtain  some  "ground  truth"  regarding  operations  in  the  quarries  from  which  recordings  of 
seismic  waves  are  available  and  verify  our  theoretical  expectations.  Such  information  is  now  being 
collected  for  the  Scandinavian  area  (S.  Mykkeltveit,  personal  communication)  and  should  be 
utilized  in  future  studies  of  this  type. 

Our  work  during  the  last  two  years  demonstrated  that  the  expected  properties  described 
above  are  indeed  exhibited  by  real  data.  In  particular,  the  factorability  concept,  a  key  feature  for 
the  kinds  of  analyses  we  propose,  is  valid  for  tiie  Pn  and  Lg  phases  recorded  at  regional  distances 
at  NORESS,  as  evidenced  by  comparisons  of  some  reconstructed  waveforms  with  actual  Pn  and 
Lg  phases  in  Figure  2.  We  have  also  found  that  Pn  waveforms  of  some  event  pairs  at  the  Titania 
mine  can  be  transformed  into  each  other  simultaneously  at  all  sensors  with  simple  common  short 
FIR  filters  (Figure  3).  For  the  same  event  pair  a  single  short  filter  designed  from  the  whole 
seismogram  increased  the  waveform  similarity,  expressed  in  tenns  of  the  site-averaged  correlation 
coefficient  from  0.04  to  near  0.7,  the  whole  set  of  arrivals.  For  other  event  pairs,  this  could  not  be 
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FIGURE  15:  Site-averaged  coherence  results  from  the  joint  factoriMtion  of  four 
E9  events  with  reconstructions  of  one  event  being  much  less  similar  to  the 
original  data  than  the  others. 
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FIGURE  16: 


Envelopes  of  band-pass  filtered  traces  Type  I  of  E9  Event. 
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FIGURE  17:  Envelopes  of  band-pass  filtered  traces  Type  II  of  Event  E9.  Even 
though  these  events  seem  to  be  noisier  than  those  if  Figure  16,  the  S/N  ratio  for 
Lg  for  these  events  is  about  the  same  as  for  E9. 
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FIGURE  18:  Coherences  demonstrating  the  high  factorization  efficiency  -  90  to 
95%  of  the  energy  accounted  for,  across  the  3C  subarray  of  NORSAR,  of  the  Pn 
phases  from  a  group  of  closely  located  FENNOLORA  profiling  shots.  What  is 
noteworthy  about  this  is  the  relatively  large  spacing,  about  5  km,  of  the  subarray 
sensors. 
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done  despite  the  fact  that  all  would  factor  well  for  individual  arrivals.  Moreover,  we  have  found 
that  events  at  the  Estonian  mine,  designated  as  E9,  fell  into  two  catagories  with  respect  to  the 
effectiveness  of  factorability.  The  Lg  phases  from  some  events,  while  still  located  at  the  same 
mine,  did  not  reconstruct  well  when  factored  with  a  group  of  others  (Figure  4).  On  closer 
examination, we  have  found  that  these  had  slightly  different  Lg  envelope  shapes,  a  feature  that  was 
very  subtle  and  could  have  gone  unnoticed.  This  may  point  to  a  different  model  composition  for 
their  Lg  phases  (and  source  mechanism). 

Summary  of  the  Processing  Results  for  Regional  Array  Data 

In  this  study  of  regional  seismic  waveforms,  extending  the  concept  of  spectral  factorability 
of  teleseismic  body  wave  data  (Filson  and  Frasier,  1972),  we  have  demonstrated  that  spectra  of 
regional  phases  can  also  be  decomposed  into  source  and  recording  site  factors,  provided  that  the 
signals  originated  from  a  limited  source  region.  Groups  of  events  located  close  to  each  other  can 
be  identified  by  the  fact  that  their  waveforms  can  be  reconstructed  from  their  spectral  source  and 
site  factors.  This  opens  the  way  for  relative  source  studies,  even  though  the  Green's  functions 
may  not  be  known  well  enough  to  model  their  waveforms  deterministically. 

In  addition  to  the  validation  of  the  factorability  concept,  it  was  also  found  that  events  within 
such  factorable  groups  could  be  further  subdivided,  according  to  relative  inter-event  coherence 
among  them.  The  events  which  show  high  coherence  probably  are  both  close  to  each  other  and 
have  very  similar  source  mechanisms,  although  their  source  time  functions  may  be  different 
enough  to  make  their  waveforms  dissimilar.  The  opposite  must  be  true  for  events  that  do  not  show 
high  coherence.  Cross-event  equalization-filtering  between  coherent  events,  using  short  time 
domain  filters  resulted  in  increases  in  waveform  similarity,  and  no  such  increase  seemed  possible 
between  event  pairs  with  low  inter-event  coherence. 
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Site-transfer  functions  between  sensors  located  at  the  extreme  ends  of  NORESS  appear  to 
be  complex,  describable  only  with  transfer  functions  with  impulse  responses  longer  than  a  second. 
Assuming  the  same  degree  of  crustal  heterogeneity  for  a  source  region,  this  may  imply  changes  in 
waveforms  similar  in  nature  over  small  displacements  of  source  location.  Further  work  needs  to 
be  done  in  exploring  the  nature  of  such  transfer  functions.  The  work  presented  demonstrates  that 
there  is  a  considerable  amount  of  relative  source  and  path  information  that  may  be  extracted  from 
regional  data,  applying  fairly  standard,  simple  coherent  time-series  processing  methods  without 
requiring  the  knowledge  or  sophisticated  modeling  of  propagation  characteristics. 
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2.0  FURTHER  ANALYSES  OF  DECONVOLUTION  RESULTS 


General  Discussion 

Waveforms  of  teleseismic  P  waves  carry  a  considerable  amount  of  information  about  the 
source  and  path,  but  such  information  is  often  hidden  by  extraneous  effects,  such  as  background 
noise,  near  receiver-scattering  (Gupta,  1989)  and  waveform  distortion,  and  near  source- scattering 
of  surface  waves  into  P  waves  (Lay,  1985).  Since  waveforms  are  complex,  it  is  not  possible  to 
see,  by  mere  visual  inspection  of  waveforms,  any  consistent,  repeating  frequency  domain  patterns 
diagnostic  of  sources  and  common  sites.  Such  patterns  exist,  nevertheless,  and  this  is  evidenced 
by  the  factorability  of  spectra  of  closely  spaced  events  observed  at  common  sets  of  sensors.  We 
are  interested  in  the  waveforms  of  P  waves  from  explosions  in  particular,  although  waveforms 
from  earthquakes  are  also  important  for  effective  discrimination. 

In  yield  estimation,  we  are  interested  in  the  particulars  of  teleseismic  P  waveforms  and  their 
effects  on  the  magnitude  estimates.  We  would  like  to  know  whether  surface  reflections  are  present 
and  what  effects  they  have  on  the  P  wave  amplitudes.  The  mere  notion  of  a  P  and  pP  "pulse"  with 
a  small,  a  few  tenths  of  a  second,  time  delay  between  them  carries  the  implication  that  the  signal 
has  some  effective  bandwidth  Be  (see  a  definition  of  this  concept  in  the  introduction  of  this  report), 
since  no  pulses  can  be  constructed  from  narrowband  data,  and  no  time  delays  of  that  order  can  be 
resolved  with  severely  band  limited  data.  On  the  other  hand,  seismic  signals  often  have  strongly 
peaked  spectra  and  small  Bg  predetermined  by  the  instrument  and  Q  along  the  paths.  Therefore, 
some  data  processing  is  required  to  recover  information  over  the  maximum  available  bandwidth 
with  acceptable  S/N  ratio  as  outlined  by  Equation  (1).  The  most  common  way  to  accomplish  this 
in  geophysics  is  to  apply  some  kind  of  signal  deconvolution  (Bakun  and  Johnson,  1973;  Douglas, 
1981;Deretal,  1987). 
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Various  deterministic  models  of  P  waveforms  from  nuclear  explosions  have  been  proposed 
and  used  in  attempts  to  resolve  multiple  arrivals.  Curiously,  'ew  attempts  were  made  to  verify  the 
consistency  of  these  models  with  the  data,  although  such  verification  is  quite  possible  using  all  the 
information  available  in  the  data.  There  are  three  basic  models  of  the  P  waves  from  explosions 
that  we  can  test: 

a)  The  most  restrictive  (P+pP)  model  assumes  that  the  P  waves  consist  of 
convolutions  of  a  sequence,  consisting  of  a  positive  P  and  a  delayed 
negative  pP  pulse  (of  identical  shapes),  convolved  with  some  time  domain 
impulse  response  functions  of  the  individual  sites,  attenuation  operator  and 
instrument  response.  Other  secondary  arrivals  are  assumed  to  be  negligible 
at  least  in  the  first  five  seconds  of  the  P  wave. 

b)  This  model  is  a  more  general  one  and  it  assumes  that  the  site  function 
concept  is  still  valid,  but  the  source  function  is  of  a  more  general  shape, 
i.e.,  it  may  contain  non-linear  distortions  of  the  P  and  pP  pulses  and 
additional  arrivals  such  as  spall.  The  spectra  Fy  (co)  of  a  suite  of  N  events 
(indexed  by  i)  recorded  at  a  set  of  M  sites  (indexed  by  /),  can  be  factored  in 
the  standard  manner. 


Fi.((y)  =  5;  iQ})Rjico) 

where  the  Si  and  Rj  are  the  source  and  recording  site  spectral  factors, 
respectively.  Generally,  it  is  advantageous  if  the  data  set  has  a  high 
redundancy  MN-M-N  and  thus,  the  number  of  total  data  far  exceeds  the 
sum  of  factors  to  be  estimated. 

Site  effects  have  been  found  to  be  quite  variable  at  some  arrays  (NORSAR, 
LASA),  thus  making  the  extraction  of  meaningful  source  information  from 
waveforms  recorded  at  individual  sensors  difficult.  At  other  array  locations 
(such  as  the  UK  arrays),  site  effects  were  found  to  be  less  important, 
apparently  due  to  the  care  in  the  selection,  using  pre-siting  surveys,  of  these 
array  sites.  It  is  reasonable  to  assume  that  the  effects  of  the  site  impulse 
response  function  may  be  reduced  by  averaging  the  data  from  numerous 
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sites,  especially  for  larger  teleseismic  arrays.  This  is  the  assumption  made 
in  the  multi-channel  deconvolution  method  developed  by  Der  et  al  (1987) 
and  the  phaseless  seismogram  methodology  of  the  UK  group  (Douglas, 
1981;  Stewart  and  Douglas,  1983).  It  has  been  shown,  on  the  other  nand, 
that  crustal  P  wave  responses  may  contain  various  converted  phases  that 
may  not  average  out  (Owens  et  al,  1987).  Nevertheless,  such  phases  are 
commonly  late  and  have  small  amplitudes  on  the  vertical  components  and 
thus,  should  not  unduly  interfere  with  the  assessments  of  early  arriving 
large  arrivals  such  as  pP. 

c)  The  third  model  is  the  most  general  one  and  does  not  even  presuppose  the 
existence  of  factorability.  This  situation  may  occur  if  we  analyze  data  for  a 
suite  of  explosions  recorded  at  a  network  of  stations  surrounding  a  test  site 
when  the  waveforms  radiated  in  the  various  directions  are  quite  different  for 
at  least  some  of  the  events  in  a  group  being  factorized.  In  this  case,  the 
spectral  factorization  of  a  suite  of  events  as  defined  above  will  not  work 
unless  all  the  events  have  similar  radiation  symmetries. 


There  are  some  statistically  easily  testable  consequences  to  each  model.  If  we  adopt  Model 
a)  by  inter-correlating  waveforms  of  pairs  of  events  using  some  best  estimates  of  P  and  pP 
parameters  (this  should  make  the  waveforms  of  a  pair  of  events  nearly  identical  at  the  same 
stations),  the  validity  of  this  signal  model  could  be  easily  tested  quantitatively.  Assuming  that  each 
seismic  trace  consists  of  a  convolution  of  a  site  impulse  response  function,  a  source  pulse  s(t)  and 
a  double  pulse  sequence  due  to  a  direct  P  and  a  pP,  P(t),  then  two  seismograms  recorded  at  the 
same  site  with  the  site  function  r(t)  for  two  different  nuclear  explosions  can  be  written  as 

tiit)=s^{t)*P^{t)*r  ii)  (10) 
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Where  *  denotes  convolution.  If  we  convolve  each  of  these  seismograms  with  the  source  pulse 
and  the  P+pP  pulse  of  the  other  event  in  the  chosen  pair  of  explosions,  then  after  this  "inter¬ 
correlation,"  the  two  events  should  produce  identical  waveforms  at  the  same  sites  (Mellman  and 
Kaufman,  1981;  Lay,  1985).  Thus, 

tj(r )  *  )  *  P^it )  =  )  *  sft )  *  Pj(t )  =  s,(r )  *S2(t )  *  P^(t )  *  P2(i  )*r(t) 

(11) 

The  similarity  of  the  two  suites  of  traces  then  could  be  tested  quantitatively  by  computing  a 
coherence  measure  such  as  that  given  in  Equation  (4)  in  the  beginning  of  this  report. 

These  coherences  should  be  close  to  unity  for  high  S/N  events  if  Model  a)  is  valid  and  the 
pP  parameters,  delays  and  relative  amplitudes  were  correct.  If  we  find  that  we  cannot  account  for 
most  of  the  power  in  the  P  waves  over  a  sufficiently  broad  frequency  band  by  this  model  using  any 
plausible  P  and  pP  parameter  combinations,  i.e.,  the  coherences  turn  out  to  be  considerably  and 
significantly  below  unity,  then  we  must  reject  Model  a).  Model  b)  requires  only  that  the  spectra 
factor  into  general  source  and  site  factors  and  that  the  reconstructions  of  the  waveforms  from  these 
factors  must  be  coherent  with  the  original  data.  For  testing  this  model  over  a  large  matrix  of 
seismograms,  we  use  the  form  given  in  Equation  (4)  in  the  introduction  of  this  report  and  we  have 
to  accumulate  the  smoothed  cross-products  of  Fourier  transforms  over  all  data  vs.  reconstruction 
pairs  for  all  the  jointly-factored  events  in  the  same  manner  (Figure  19a).  In  measuring  this 
coherence,  we  must  also  adjust  the  degrees  of  freedom  so  that  we  take  into  account  the  fact  that  we 
have  already  estimated  the  site  and  source  factors  from  the  same  data.  It  is  also  possible  to  test  the 
data-reconstruction  similarity  only  over  a  single  event  for  all  sites  or  a  single  site  for  all  available 
events  by  performing  the  additions  only  over  sites  or  events,  respectively,  as  it  was  done  in 
Section  1.0  for  regional  events.  Failure  to  get  high  values  of  C  in  this  test  would  eliminate  Model 
b)  also  for  the  given  set  of  events.  As  opposed  to  simple  trace  comparisons,  the  coherence 
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FIGURE  19:  Methods  for  testing  the  efficiency  (portion  of  energy  accounted 
for),  a)  By  spectral  factorization,  b)  By  the  P+pP  model. 
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measurement  described  above  is  not  as  sensitive  to  minor  differences  in  the  assumed  source 
models  such  as  those  most  widely  used  (von  Seggem-Blandford,  1972;  Mueller-Murphy,  1971; 
and  Lay,  1985).  Thus,  any  results  would  not  be  affected  if  one  did  not  use  the  same  source  pulse 
models  as  the  authors  of  the  original  papers  being  evaluated  below  and  this  approach  allows  us 
more  flexibility. 

In  case  Model  a)  is  found  to  be  not  valid,  we  can  no  longer  assume  that  the  spectral  (or 
cepstral)  methods  will  yield  reliable  or  meaningful  estimates  of  "pP"  parameters,  since  multiple 
arrivals  consisting  of  P,  pP,  spall  and  other  arrivals  will  affect  the  spectral  and  cepstral  analyses  in 
very  complex  and  generally  intractable  ways,  even  though  pP  phases  may  be  consistently  present 
at  all  stations.  Moreover,  shapes  of  absolute  values  of  spectra  discarded  all  the  phase  information 
in  the  signals  already  and  are  not  uniquely  invertible.  For  instance,  the  spectra  are  unchanged  by 
reversing  the  polarity  or  time  sequence  of  the  time-series  or  filtering  them  with  dispersive  filters 
with  flat  amplitude  response  that  shift  the  phase  only.  In  the  case  the  P  wave  spectra  do  not  factor 
over  all  the  sensors,  (case  c),  it  is  hard  to  talk  about  any  meaningful  "pP"  consistent  over  the 
network  of  stations.  If  Model  b'  ’S  valid,  then  we  can  still  estimate  the  source  time  function.  In 
the  case  of  c)  not  much  can  be  done  except  that  we  may  still  use  small  arrays  to  stack 
decovolutions  at  the  individual  sensors  over  which  factorability  holds  to  look  at  the  source  from 
various  directions. 

It  must  be  pointed  out  that,  although  our  understanding  of  the  physical  processes  that 
follow  nuclear  explosions  is  poor,  it  is  not  very  likely  that  the  P+pP  model  could  adequately 
characterize  P  waves  from  nuclear  explosions.  Non-linear  processes,  such  as  spalling,  strain 
release,  block  motion  and  screening  of  P  waves  by  zones  of  damaged  rock  (McLaughlin  et  al, 
1988;  Barker  et  al,  1990;  Douglas  and  Hudson,  1990)  are  likely  to  result  in  complex  waveforms 
resulting  from  multiple  secondary  arrivals  besides  the  pulse  that  can  be  interpreted  as  pP  (Douglas, 
1981).  Despite  this,  it  is  still  likely  that  the  sizes  of  first  P  pulses,  if  they  can  be  isolated,  will  be 
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relatively  uncontaminated  and  can  be  used  for  measuring  yields  (Douglas  et  al  1987;  Stewart  and 
Douglas,  1983).  Alternatively,  corrections  to  the  conventional  mt),  based  on  the  deconvolution 
results,  could  be  conceivably  derived. 

Arguments  in  Favor  of  Deconvolution  Results 


Despite  the  sometimes  counter-intuitive  and  complex  nature  of  the  waveforms  derived  from 
multi-channel  deconvolution,  physical  arguments  can  be  made  for  the  claim  that  these  represent  real 
features  of  the  source-related  radiation.  These  arguments  can  be  summarized  as  follows: 

a)  Explosions  in  hard  rock  with  flat  overlying  topography  generally  have  a 
good  "pP"  feature. 

b)  The  cratering  explosion  at  Kazakh  did  not  have  a  good  pP. 

c)  Explosions  in  areas  with  complex  topography  (Degelen  and  Algeria)  have 
poor  and/or  azimuthally  strongly  varying  pP  phases. 

d)  Combined  t*  and  pP  effects  estimated  from  spectral  Q  studies  and 
deconvolutions  result  in  magnitude-yield  relationships  that  seem  to  agree 
with  those  observed  for  explosions  with  known  yield  (see  detailed 
discussion  later  in  this  report). 

e)  Small,  middle-sized  earthquakes  deconvolved  in  the  same  manner  (or 
equivalently  broadband  seismograms)  often  show  clear,  simple  pP  and  sP 
pulses  one  would  expect  for  point  sources.  This  shows  that  the 
deconvolution  method  does  not  give  spurious  waveforms  for  simple 
sources. 

0  The  spectral  factorization  model  accounts  for  90  to  95%  of  the  signal  energy 
over  broad  frequency  bands  (1-6  Hz  where  the  S/N  allows  it)  in  the  multi¬ 
sensor  seismogram  assemble  why  source  estimation  procedures  based  other 
models  account  for  much  less. 
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Figure  20  illustrates  some  typical  deconvolution  results  from  a  few  major  test  sites  referred 
to  above  and  Figure  21  illustrates  the  validity  of  item  f),  above. 

Tests  of  Some  pP  Results  from  the  Literature 

Estimates  for  pP  parameters  for  a  set  of  Pahute  Mesa  explosions  were  computed  by  Lay 
(1985)  and  Murphy  (1989).  The  methods  applied  were  different:  Lay  used  a  version  of  the  inter¬ 
correlation  method  developed  by  Mellman  and  Kaufman  (1980)  which  maximizes  the  time  domain 
correlation  coefficient  between  suites  of  inter-correlated  traces,  while  Murphy  used  the  shapes  of 
network-averaged  spectra  to  estimate  the  pP  parameters.  In  Murphy's  work  the  spectral  estimates 
are  computed  in  a  theoretically  quite  intractable,  unconventional  fashion  by  using  sets  of  bandpass 
filters  and  picking  the  maximum  amplitude  within  the  first  five  seconds  of  the  filtered  signal. 
Appendix  C  describes  an  implementation  of  the  MCD  method  in  the  X  winds  environment.  We 
see  no  advantages  for  such  unproven  procedures  in  view  of  the  many  existing  spectral  estimation 
methods  with  well-understood  characteristics  (Marple,  1987).  Moreover,  the  physical  meaning  of 
such  a  spectral  measure  is  not  clear  and  is  further  obscured  by  elaborate  and  non-physical 
correction  and  station  averaging  schemes  applied,  allegedly,  to  compensate  for  path  t*,  station 
effects  and  pP.  Instead  of  attempting  to  analyze  these  physically  and  statistically  intractable 
procedures  we  shall  concentrate  on  checking  some  of  their  consequences  with  regard  to  pP 
parameters. 

In  order  to  test  the  validity  of  Murphy's  results  for  pP  parameters  for  Pahute  Mesa  events 
(Murphy,  1989)  from  network-averaged  spectra  (NAS),  we  have  used  the  inter-correlation  method 
of  Lay  (1985),  implemented  as  shown  in  Figure  19.  Since  the  two  source  pulse  models  are  not 
very  different,  our  conclusions  are  not  affected  much  by  the  fact  that  we  used  the  von  Seggem- 
Blandford  pulse. 
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FIGURE  20:  Some  typical  deconvolution  results  for  various  test  sites. 
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FIGURE  21:  Data  and  reconstructions  using  the  spectral  factorization  method. 
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We  have  attempted  to  verify  Murphy's  pP  parameters  by  inter-correlating  multi-sensor 
NORSAR  data  for  the  same  Pahute  Mesa  explosions.  Using  NORSAR  data  had  the  advantages 
that  the  same  events  were  digitally  recorded  over  the  array  and  that  the  instrumental  response  at 
NORSAR  results  in  waveforms  truly  broadband  in  nature,  as  opposed  to  the  narrowband 
instruments  of  the  WWSSN  network,  utilized  by  Murphy  (1989)  and  Lay  (1985)  where  the 
bandwidth  was  further  reduced  by  hand-digitization.  Inter-correlation  reintroduces  the  phase 
information  in  the  signals  that  was  discarded  by  the  computation  of  network-averaged  spectra  and 
provides  a  conclusive  test  for  the  implicitly  assumed  P-i-pP  model.  We  approximately  compensated 
the  source  waveforms  for  source-scaling  by  using  a  cube-root  scaled  von  Seggem-Blandford 
model,  deconvolving  the  source  waveform  appropriate  to  one  of  the  events  and  applied  to  it  the  one 
that  was  appropriate  to  the  other.  In  accordance  with  the  time  windows  used  in  both  inter¬ 
correlation  and  NAS  methodology,  we  computed  our  correlation  coefficients,  using  a  five  second 
window  encompassing  only  the  first  part  of  the  P  wavetrains  at  the  various  sensors.  Figure  22 
shows  some  of  the  waveform  inter-correlation  results  for  two  pairs  of  nuclear  explosions.  The 
wavetrains  that  resulted  from  the  procedure  for  the  pairs  of  events  as  seen  at  the  same  sensors  are 
not  similar  at  all,  not  even  in  the  first  five  seconds,  thus,  indicating  that  either  Model  a)  or  the  pP 
parameters  of  Murphy  (1989),  or  both,  are  not  correct.  Tabulating  the  time  domain  correlation 
coefficients  (Table  III)  demonstrates  the  fact  that  the  pP  parameters  derived  by  Murphy  from 
network  averaged  spectra  and/or  the  assumption  of  the  simple  P+pP  model  cannot  be  reconciled 
with  the  NORSAR  data.  Contrary  to  the  expectation  for  a  correct  model,  the  correlation 
coeuic’ents  do  not  increase  significantly  after  inter-correlation  and  in  one  case  they  decrease.  In 
computing  the  time  domain  correlation  coefficients,  we  have  allowed  for  the  possibility  that  the 
event  waveforms  may  not  have  been  lined  up  in  time  properly  and  used  a  five  second  time  window 
for  estimating  the  correlation  coefficients  in  order  to  conform  with  the  window  used  by  Murphy  for 
computing  the  network-averaged  spectra.  We  have,  thus  tried  to  find  the  best  time  delay  with  the 
highest  correlation  coefficient  by  shifting  all  the  waveforms,  using  the  same  delays  for  all  traces  for 
one  event.  The  resulting  best  correlation  coefficients  are  presented  in  this  report.  To  estimate  how 
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FIGURE  22:  Waveform  inter-correlation  results  at  NORSAR  using  Murphy’s 
(1989)  pP  parameters  for  Pahute  Mesa  events.  The  waveforms  at  the  same 
sensors  should  be  identical,  but  they  are  not  indicating  that  the  P+pP  model 
and/or  the  pP  parameters  are  invalid. 
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TABLE  III:  Intercorrelation  Results  with  pP  Parameters  Derived  from  Spectra 


Event  Pair 

Average  Correlation 

Coefficient 

Before 

After 

Muenster/Inlet 

0.234 

0.349 

Almendro/Mast 

0.570 

0.619 

Tybo/Inlet 

0.577 

0.552 

Muenster/Mast 

0.336 

0.521 
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much  power  is  accounted  for  by  the  pP  parameters  and  the  P+pP  model,  we  have  computed  a 
heavily-smoothed,  site-averaged  coherence  estimate  between  the  inter-correlated  P  waveforms  at 
NORSAR  for  a  few  pairs  of  events.  The  time  window  used  was  again  five  seconds  in  length  for 
consistency  with  both  the  NAS  and  inter-correlation  methodology.  The  result  for  three  pairs  of 
events  are  shown  with  the  dashed  lines  in  Figure  23,  in  the  middle  row.  Compared  with  the 
factoring  results,  which  account  for  more  than  95%  of  the  energy  of  P  waves,  the  results  from 
network  averaged  spectra  account  for  much  less.  For  comparison,  we  also  show  a  similar  analysis 
for  data  vs.  reconstruction  pairs  after  multi-channel  deconvolution  for  a  sample  of  Kazakh  events 
on  the  right  of  the  figure  (NTS  events  would  give  similar  results,  since  typical  correlation 
coefficients  between  data  and  reconstructions  are  near  0.9)  at  EKA,  taking  into  account  that  the 
degrees  of  freedom  in  this  case  are  decreased  by  the  fact  that  the  site  and  source  factors  were 
estimated  from  the  same  data.  We  have  set  the  degrees  of  freedom  in  both  cases  near  200,  thus 
both  are  comparable.  Any  small  differences  between  the  von  Seggern-Blandford  and  Mueller- 
Murphy  models  should  also  cancel  in  this  analysis  for  the  Muenster/Mast  pair,  since  coherence 
analyses  allow  for  the  existence  of  simple  transfer  functions.  The  comparison  shows  that  spectral 
factoring  accounts  for  about  90%  of  the  total  signal  energy  even  in  this  very  wide  (0.5  -  6)  Hz 
band  in  the  first  five  sends  of  the  P  waves.  We  must  note,  however,  that  even  if  we  had  used  the 
customary  window  of  25  seconds  for  coherence  estimation,  the  spectral  factorization  would  have 
still  accounted  for  about  the  same  percentage  of  total  energy.  These  findings  indicate  that  the  P-t-pP 
model  and/or  the  pP  parameters  derived  by  Murphy  for  Pahute  Mesa  are  not  correct.  We  believe 
that  the  cause  of  the  poor  performance  of  the  inter-correlation  and  NAS  is  due  to  the  fact  that 
source  waveforms  of  nuclear  explosions  from  Pahute  Mesa  are  more  complex  and  cannot  be 
explained  by  the  simplistic  P+pP  model. 

Lay  (1985),  using  a  WWSSN  data  set  similar  to  that  of  Murphy  (1989),  also  estimated  a 
set  of  pP  parameters  for  essentially  the  same  set  of  Pahute  Mesa  explosions  using  the  inter¬ 
correlation  method.  In  addition  to  the  high  t*  along  the  paths  from  NTS  (Der  et  al,  1985)  and  the 
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FIGURE  23:  Portion  of  P  wave  energy  in  Pahute  Mesa  events  accounted  for  at  NORSAR  by  Lay’s  (top  row), 
Murphy’s  pP  parameters  (second  row)  all  shown  by  dashed  lines.  The  energy  accounted  for  in  Pahute  Mesa 
events  by  spectral  factorization  (multi-channel  deconvolution)  is  shown  by  solid  lines  in  the  first  two  rows.  The 
bottom  graph  shows  the  average  value  of  the  same  for  Kazakh  events  obtained  by  spectral  factorization. 


narrowband  response  of  the  WWSSN  instruments,  the  effective  bandwidth  utilized  in  his  analysis 
was  further  reduced  by  reapplying  a  source  spectrum  with  narrowband  characteristics  (Hartzell  et 
al,  1983).  The  inter-correlations  resulted  in  some  extremely  narrowband  wavelets  which  are 
essentially  strongly-tapered  sinusoidal  transients.  Typically,  he  obtained  time  domain  correlation 
coefficients  above  0.8,  but  only  a  few  exceeded  0.9.  The  waveform  examples  in  the  1985  paper 
are  nonrepresentative  since  mostly  those  with  the  highest  coefficients  are  shown.  These  relatively 
high  coefficients,  on  the  other  hand,  are  to  be  expected  for  even  random  wavelets  with  narrow 
bandwidth,  since  sample  correlation  coefficients  for  severely  band-limited  data  (low  degrees  of 
freedom)  are  strongly  biased  towards  high  values  (Der  et  al,  1983;  Bendat  and  Piersol,  1966). 

Moreover,  some  waveforms  have  been  shifted  individually  in  time  to  achieve  the  "best  fit" 
in  such  a  manner  that  the  first  arrivals  no  longer  coincide.  We  are  unable  to  think  of  any 
justification  for  such  a  procedure,  since  it  conflicts  with  the  basic  signal  model  assumed.  The 
optimization  in  Lay's  study  was  performed  by  fixing  the  parameters  of  one  event  and  determining 
the  pP  of  the  other  relative  to  it.  Lay's  results  were  non-unique  because  they  depended  strongly  on 
the  choices  of  the  pP  parameters  for  the  "master"  event.  Thus,  there  are  several  sets  of  mutually 
incompatible  results  for  the  same  events. 

We  have  used  the  same  NORSAR  data  for  testing  the  best  pP  parameters  (those  with  the 
highest  correlation  coefficients)  derived  by  Lay.  The  results  (top  row  in  Figure  4)  again  show  that 
the  P+pP  model  does  not  account  for  a  significant  portion  of  the  total  P  wave  energy  in  the  first 
five  seconds  at  the  NORSAR  array  although  it  does  generally  better  than  Murphy's  parameters. 
This  is  not  surprising  since  Lay's  methodology  takes  the  signal  phase  information  into  account 
while  Murphy's  does  not  and  was  designed  to  maximize  the  correlation  coefficient.  The  coherence 
analyses  only  used  the  first  five  seconds  of  the  signal  to  conform  with  the  quite  reasonable 
assumption,  made  by  both  Murphy  and  Lay,  that  the  latter  part  of  the  signal  mostly  contains 
scattered  waves  and  has  nothing  to  do  with  the  source  processes.  From  the  preceding,  we  must 
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conclude  that  we  find  no  support  from  NORSAR  data  analyses  for  the  P+pP  signal  model  or  the 
pP  parameters  derived  from  it. 

In  order  to  test  the  pP  results  of  Murphy  and  Lay  over  a  larger  network  addition  to 
NORSAR  data,  we  also  acquired  sets  of  hand-digitized  WWSSN  data,  which  contain  the  same 
seismograms  as  those  used  by  Murphy  (1989)  and  Lay  (1985).  Performing  inter-correlations  on 
suites  of  waveforms  using  the  two  sets  of  pP  parameter  estimates,  we  have  arrived  at  the  results 
shown  in  Table  IV.  The  inter-correlations  generally  increased  the  correlation  coefficients  for  both 
sets  of  pP  parameters.  In  order  to  test  the  significance  of  these  increases,  we  have  computed  the 
effective  mean-square  bandwidth  defined  in  Equation  (2)  for  our  data  sets.  We  estimated  the 
confidence  limits  for  the  correlation  coefficients  using  the  computed  values  of  Be  and  the 
algorithms  in  Bendat  and  Piersol  (1966).  Only  in  one  case  is  the  increase  in  correlation  coefficient 
is  significant  statistically,  in  two  cases  it  is  marginal,  close  to  the  95%  confidence  limit,  and  in  four 
cases  it  is  not  significant.  It  is  interesting  to  note  that  the  effective  bandwidth  for  such  data  is  only 
0.3  Hz,  and  after  we  reimposed  the  source  frequency  spectrum  the  second  time  in  the  manner  of 
the  original  inter-correlation  paper  (Lay,  1985)  this  bandwidth  further  decreased  to  0.15-0.18  Hz. 
The  small  Be  are  the  results  of  the  low  path  Q  for  teleseismic  P  waves  from  NTS  (Der  et  al,  1985) 
and  the  narrow  spectral  response  of  short  period  WWSSN  instruments.  The  above  results  should 
vividly  illustrate  the  dangers  of  applying  time  domain  analysis  methods  to  severely  band-limited 
data. 
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TABLE  IV:  Intercorrelation  Results  with  WWSSN  Data  Using 
Published  pP  Parameters  for  Pahute  Events 


Event  Pair 

Method 

Corr.  Coeff. 
before 

Corr.  coeff. 
after 

Significant? 

Estuary/Mast 

NAS 

0.378 

0.458 

No 

OIC 

0.378 

0.588 

No 

OICS 

0.398 

0.697 

Marginal 

Mast/Tybo 

NAS 

0.639 

0.660 

No 

OIC 

0.639 

0.660 

No 

OICS 

0.659 

0.681 

No 

Estuary/Tybo 

NAS 

0.563 

0.670 

Marginal 

OIC 

0.563 

0.730 

Yes 
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The  results  discussed  above  prove  that  the  P+pP  model  does  not  explain  a  significant 
portion  of  the  signal  energy  and  that  the  sometimes  unexplained  complexities  in  the  source 
waveforms  shown  in  the  deconvolved  P  waves  are  needed  for  explaining  the  data  and  accounting 
for  the  energy  in  the  signal  in  a  broader  frequency  band.  It  appears  that  the  pP  parameters  derived 
by  Murphy  (1989)  and  Lay  (1985)  cannot  be  substantiated  and  are  probably  artifacts.  Any  results 
from  spectral  estimation  for  pP  parameters  can  also  be  questioned  on  the  basis  that  spectra  are  not 
invertible  unless  we  know  a  priori  that  the  P+pP  model  is  valid.  The  results  of  Lay  (1985)  fit  the 
data  better,  but  the  large  residual  power  unexplained  by  the  P+pP  model  from  both  studies 
indicates  that  P  waves  from  Pahute  Mesa  are  quite  complex  and  that  the  complex  source  time 
functions  resulting  from  the  deconvolutions  are  necessary  to  account  for  the  broadband  spectral 
structure  of  the  data.  The  bandwidth  reduction  associated  with  Lay's  (1985)  procedures  make 
them  fail  most  tests  for  statistical  significance.  Therefore,  for  Pahute  Mesa  explosions  the  P+pP 
model  fails  the  test.  The  P+pP  model,  on  the  other  hand,  may  be  reconciled  with  the  Kazakh  data 
better,  since  the  deconvolved  source  functions  of  these  events  are  simpler,  although  possible 
additional  smaller  "spall"  phase  and  other  arrivals  may  also  be  present  in  many  of  the 
deconvolutions  besides  the  "pP".  Nevertheless,  since  the  deconvoh  tion  results  must  be  inspected 
first  to  confirm  these  notions,  we  see  no  reason  to  use  the  other  two  methods,  we  have  just 
evaluated  for  Kazakh  data  either. 

Deconvolution  and  its  Consequences  with  Respect  to  Magnitude-Yield 

All  the  discussions  about  the  nature  of  P  waves  is  of  academic  interest  only  unless  the 
secondary  arrivals  have  an  appreciable  effect  on  the  wave  amplitudes  and  thus  the  yield  estimates. 
The  motivation  for  developing  our  multi-channel  deconvolution  method  was  to  gain  information 
about  seismic  sources,  most  notably,  nuclear  explosions  and  to  improve  yield  estimation 
procedures.  In  the  past,  our  ideas  about  the  waveforms  radiated  from  nuclear  explosions  were 
based  more  on  theoretical  preconceptions  than  on  observational  facts.  Therefore,  we  still  do  not 
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know  what  the  radiated  waveforms  from  nuclear  explosions  would  look  like  if  we  had  senscrs 
embedded  in  the  medium  outside  the  elastic  radius.  The  results  from  our  deconvolutions  revealed 
that  not  all  nuclear  explosions  seem  to  have  a  "pP"  reflection.  Explosions  at  Shagan,  Sinkiang  and 
Novaya  Zemlya  generally  have  a  pronounced  pP,  while  those  at  Pahute  Mesa,  Ahaggar  and 
Degelen  do  not.  All  work  on  "pP"  is  only  of  academic  interest  unless  the  result  can  be  used  for 
yield  estimation.  The  question  is  whether  we  still  could  assume  that  the  first  pulse  in  the 
deconvolved  wavetrain  still  scales  with  yield,  such  that  the  amplitude  of  the  total  P  wave  is 
approximately  determined  by  the  superposition  of  the  wavelets  from  this  first  pulse  and  a  pP.  If 
so,  then  the  pP  will  have  a  strong  effect  on  P  wave  amplitudes  and  hence  on  yield  estimates.  We 
can  expect  then  that  for  test  sites  where  the  pP  is  generated  the  wave  amplitude  will  generally  be 
reinforced.  Thus,  in  addition  to  the  biases  due  to  Q,  we  may  have  to  deal  with  biases  associated 
with  pP. 

In  our  previous  work,  we  have  done  some  simulations  of  magnitude-yield  (using 
conventional  mb)  relationships  by  using  the  von  Seggem-Blandford  (1972)  source  model  in 
conjunction  with  a  constant  "effective"  t*  operator,  a  WWSSN  short  period  instrument  response 
and  the  assumption  that  the  first,  direct  pulse  scales  directly  with  yield.  Due  to  the  frequency 
spectra  of  the  von  Seggern-Blandford  source  the  attenuation  effect  was  strongly  yield-dependent 
resulting  in  curved  rather  than  straight  line  magnitude-log  (yield)  relationships.  In  this,  we  have 
followed  similar  work  done  previously  to  justify  the  basis  for  magnitude  yield  formulas  using 
some  physical  model.  In  the  simulations  we  have  assumed,  based  on  generalized  deconvolution 
results,  that  only  Shagan  events  had  a  pP  and  that  the  t*  differences  conformed  with  rhose  tound 
from  spectral  studies  of  P  waves  (Der  et  al,  1985).  The  waveforms  resulting  from  the  simulations 
were  subjected  to  the  usual  body  wave  calculation  procedures  and  were  plotted  in  Figure  24.  The 
two  sets  of  curves  resulted  by  assuming  that  the  recording  station  was  situated  in  a  shield  area  (a) 
or  in  a  tectonic  area  (b)  with  higher  mantle  attenuation.  The  attenuation  related  source  biases  would 
be  more  visible  with  recording  stations  situated  in  shields,  since  tnrough  a  second  pass  through 
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high  Q  mantle  the  differentials  in  the  wave  amplitudes  would  be  preserved.  In  the  case  of  stations 
over  a  low  Q  mantle,  the  differential  would  be  decreased  due  to  the  relative  removal  of  the  high 
frequency  components  of  signals.  The  curves  in  Figure  A1  did  not  have  any  confirmation  in  the 
real  data  at  the  time  they  were  produced,  in  fact  the  authors  of  the  original  report  (Der  et  al,  1987) 
were  themselves  quite  skeptical  of  the  relevance  of  such  results.  According  to  these  figures,  the 
combined  effect  of  Q  and  pP  effects  would  result  in  very  large  biases  that  increased  with 
decreasing  yields.  There  would  be  a  negative  offset  between  the  Degelen  and  Shagan  test  sites  and 
the  NTS-Shagan  bias  would  be  much  higher  than  that  anyone  would  have  thought  to  be 
reasonable.  On  the  other  hand,  no  data  to  refute  such  predictions  were  available  either,  there  were 
no  compiled  measurements  of  the  kind  for  Vv  iVSSN  stations  alone,  these  stations  were  (and  still 
are  not)  classified  into  low  and  high  Q  stations,  and  finally,  there  were  no  yields  announced  for 
most  Soviet  nuclear  explosions.  The  then  available  body  wave  magnitude  data  were  compounded 
from  measurements  on  a  variety  of  instruments  and  networks,  often  using  empirical  relative  scalar 
additive  terms  for  equalizing  the  body  wave  magnitudes  among  the  various  types  of  instruments, 
which  cannot  possib’y  take  care  of  ad  the  possible  physical  effects  of  attenuation  and  pP.  The 
simple  linear  empirical  magnitude-yield  relationships  characterized  with  the  intercept  and  slope 
parameters  and  with  constant  "bias"  and  offsets  among  test  sites  thus,  had  no  physical 
justifications  based  on  any  physical  model. 

We  have  also  repeated  the  same  kind  of  calculations  using  the  Mueller  and  Murphy  source 
model  and  a  shield  station  assumption  (Mueller  and  Murphy,  1971;  Murphy,  1977).  This  model 
differs  from  the  VS  B  model  in  that  the  source  depends  on  the  depth  explicitly  and  in  that  the  source 
time  scaling  is  less  dependent  on  yield.  Results  of  some  of  these  calculations  are  shown  in  Figure 
25.  Note  the  large  NTS-Kazakh  bias  and  the  Degelen  negative  offset  relative  to  Shagan.  In 
contrast  to  Figure  25,  the  Mg-mb  lines  have  little  curvature  and  the  biases  increase  only  slightly 
with  decreasing  yield.  The  latter  features  are  more  in  agreement  with  regression  results  on  real  data 
reported  in  the  literature. 
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Von  Seggern-BIandford  Model 


LoglO(YLD  dt) 


FIGURE  24:  Simulated  WWSSN  mb  for  the  three  test  sites  obtained  by  using  the 
source  model  of  von  Seggern  and  Blandford  (1972),  t*  derived  from  P  spectra 
and  pP  parameters  derived  from  deconvolution.  Note  the  large  NTS-Kazakh  bias 
and  the  Degelen  negative  offset  relative  to  Shagan.  The  Ms-nib  lines  have  a 
pronounced  curvature  and  the  biases  increase  with  decreasing  yield.  The  two  sets 
of  curves  are  for  stations  located  in  shield  and  tectonic  areas  respectively. 
Identifical  source  media  were  assumed  for  simplicity. 
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Recently,  Jih  (1990)  analyzed  a  large  magnitude  data  set  for  WWSSN  stations  and  some 
other  station  recordings  that  were  transformed  to  have  the  WWSSN  response.  He  performed 
regressions  on  announced  U.S.  yields  and  the  recently  announced  Soviet  yields  listed  by  Vergino 
(1989).  In  these  regressions,  the  censoring  effects  on  yields  and  the  magnitudes  were  fully 
considered.  The  resulting  magnitude-yield  regression  lines  are  shown  in  Figure  26.  The 
similarities  between  this  figure  and  Figures  24  and  25  can  be  noted.  The  NTS-Shagan  bias  is 
large,  larger  than  one  could  explain  with  t*  differentials  estimated  from  spectra  (Der  et  al,  1985).  It 
is  also  magnitude  dependent  and  there  is  a  Degelen-Shagan  offset.  These  general  features  agree 
with  the  predictions  in  Figures  24  and  25.  Before  jumping  to  conclusions,  we  must  add  a  note  of 
caution.  The  slopes  in  Figures  24  and  25  are  not  the  same  as  those  obtained  by  Jih  (1990)  from 
the  data,  and  although  the  offsets  off  the  three  curves  are  correct  in  sign,  the  differences  are  not  the 
same.  Some  of  these  discrepancies  may  be  appropriated  to  the  fact  that  the  assumptions  and  the 
source  models  were  highly  idealized  in  the  simulations.  For  instance,  actual  deconvolutions  of 
some  Degelen  and  NTS  explosions  do  show  pP-like  downswings  and  some  Shagan  explosions  do 
not  have  pronounced  "pP"  arrivals.  This  could  make  the  separation  of  these  populations  less 
pronounced.  Besides,  real  data  contains  other  arrivals  as  well. 

Jih's  results  must  also  be  further  evaluated  to  reduce  the  possibility  of  network  bias 
contamination  of  the  conclusions.  To  do  this,  we  must  strive  for  a  near-uniform  distribution  of 
stations  on  the  focal  sphere  and  the  uniformity  of  the  station  responses.  Jih  used  data  from  other 
systems  as  well  and  he  transformed  the  traces  to  WWSSN  short  period  instrument,  prior  to  making 
readings  for  magnitude  calculations.  With  these  added  data,  Jih  already  obtained  a  fairly  well- 
distributed  network  with  a  uniform  instrument  response. 

We  must  be  careful  not  to  draw  hasty  conclusions  from  Jih's  results  and  the  similarities 
between  our  earlier  simulation  results  and  those  of  Jih  (1990)  do  not  necessarily  prove  that  either 
of  them  are  right,  or  that  we  fully  understand  the  physical  factors  underlying  magnitude-yield 
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FIGURE  25:  Simulated  WWSSN  mb  for  the  three  test  sites  obtained  by  using  the 
source  model  of  Mueller  and  Murphy  (1971),  t*  dereived  from  P  spectra  and  pP 
parameters  derived  from  deconvolution.  Note  the  large  NTS*Kazakh  bias  and  the 
Degelen  Negative  offset  relative  to  Shagan.  The  Ms-mb  lines  have  little  curvature 
and  the  biases  increase  on'y  slightly  with  decreasing  yield.  The  latter  features  are 
more  in  agreement  with  the  data.  Identical  source  media  were  assumed  for 
simplicity. 
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FIGURE  26:  Jih’s  regression  results  for  NTS  and  Kazakh  (Degelen  and  Shagan) 
events  using  maximum  likelihood  analyses  of  actual  and  simulated  P  waves 
(recorded  at  other  stations  transformed  to  WWSSN  response).  Note  the  offset  of 
Degelen  relative  to  Shagan  and  the  large  NTS-Shagan  difference  that  increase 
with  decreasing  mb.  Note  the  similarities  of  these  lines  to  the  general  features  in 
Figure  B1  and  especially  to  those  in  Figure  25. 
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relationships.  Nevertneless,  the  simulation  results  point  to  the  dangers  of  the  narrow  empiricism 
that  characterized  most  of  the  studies  of  magnitude-yield  relationships,  and  yield  estimation 
procedures.  Clearly,  we  cannot  be  confident  of  the  validity  of  our  yield  estimation  procedures  until 
we  close  the  gap  between  the  simplistic  linear  magnitude-yield  formulas  differing  only  in  their 
"intercepts"  and  the  physical  models  that  predict  yield-dependent  bias,  effects  of  mantle  Q  under 
each  station  and  a  pP  effect.  Both  cannot  be  right  and  have  to  be  reconciled  by  modifying  the 
physical  models  and  by  defining  magnitude  measures  that  depend  on  the  physical  factors  in  a  more 
consistent,  meaningful  way.  The  standard  body  wave  magnitude  does  not  qualify  (Der  et  al,  1981) 
since  the  interactions  among  wave  "periods,"  the  instrument  response  and  divisions  by  period  tend 
to  obscure  the  physical  processes,  attenuation,  source  scaling  and  pP,  in  such  a  way  that 
magnitudes  will  not  depend  on  any  of  these  factors  in  any  clear,  invertible  fashion.  For  instance,  it 
can  be  shown  by  simulations  that  in  some  situations  strong  attenuation  of  P  phases  may  even 
increase  the  body  wave  magnitude  (Der  et  al,  1981).  Using  other  measures  of  body  wave 
magnitudes  such  as  true  spectral  mb-s  may  be  less  affected  by  pP  than  implied  above.  On  the  other 
hand,  band-pass  filtered  P  wave  maximum  amplitudes  used  in  computing  the  network-averaged 
magnitudes  (Murphy,  1989)  may  still  behave  similarly  to  the  conventional  mb,  although  some  of 
the  problems  with  period  "correction"  and  division  by  the  period  are  avoided,  and  thus,  such 
measures  are  still  vastly  preferable,  on  physical  grounds,  to  conventional  mt,. 
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3.0  STUDIES  OF  SPATIO-TEMPORAL  RESOLUTION  OF  SOURCE 
INVERSIONS 

(Jeneral  Discussion 

The  time  domain  method  for  the  derivation  of  source  and  path  parameters  from  teleseismic 
recordings,  direct  modeling  or  least  squares  inversion  of  the  observed  seismic  waveforms  is 
popular  and  widely  used.  The  seismological  literature  contains  hundreds  of  papers  of  this  kind, 
i’lie  objective  of  the  method  is  to  obtain  optimum  fits  between  the  observed  waveforms  and  the 
ynthetic  seismograms;  the  quality  of  fits  is  often  Judged  by  visual  inspection  or  by  optimizing 
some  objective  function,  most  commonly  minimizing  the  RMS  differences  between  the  data  and 
synthetic  seismograms  (Nabelek,  1984;  Wallace  et  al,  1981).  Some  studies  utilize  elaborate  least- 
squtu'es  inversion  schemes  with  complex  constraints  for  estimating  the  rupture  history.  Such 
methods  (applied  mostly  to  data  sets  that  also  include  strong  motion  data),  were  described  for 
example  by  Hartzell  and  Heaton  (1983),  Hartzell  (1989),  and  Hartzell  and  lida  (1990).  At  the 
point  of  optimum  time  domain  fit  it  is  assumed  that  the  parameters  used  in  computing  the  synthetic 
seismograms  are  the  best  obtainable  estimates  for  those  of  the  source  and  path.  The  data  utilized 
are  diverse;  they  may  include  long  period  teleseismic  waveforms,  teleseismic  short  period 
waveforms  and  strong  motion  data.  Although  a  narrative  of  the  procedures  leading  to  this  best 
solution  is  commonly  presented,  this  is  usually  not  followed  by  assessments  of  the  possible 
uncertainties  and  the  limits  of  spatio-temporal  resolution  by  taking  into  consideration  the  relative 
sizes  of  residual  errors.  The  commonly  missing  resolution  or  error  estimates  should  be  part  of  any 
study.  In  this  report,  we  shall  discuss  mostly  the  inversion  of  teleseismic  long  period  body 
waveform  data,  but  the  general  principles  should  apply  to  all  inversion  work. 

Upon  inspection  of  the  results  of  most  waveform  inversion  studies,  it  is  obvious  that  the 
synthetic  and  real  seismograms  do  not  match  perfectly  and  that  the  differences  are  mostly  in  the 
low  amp;:;  jdc,  hign-frequency  details  of  the  waveforms.  The  synthetic  seismograms  generally 
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also  contain  less  high-frequency  energy  than  the  original  data  even  if  we  allow  for  the  presence  of 
background  noise.  The  inverted  source  time  functions  presented  often  contain  details  of  a  few 
seconds  duration  that  depend,  for  their  adequate  characterization,  on  the  same  high-frequency 
signal  components  that  were  mismatched  in  the  waveform  inversion  or  modeling  process.  The 
time  durations  of  the  body  waves  observed  at  teleseismic  stations  often  greatly  exceed  those 
claimed  for  the  source,  often  there  is  a  somewhat  arbitrary  cutoff  in  the  length  of  time  windows 
beyond  which  no  attempt  is  being  made  to  fit  the  synthetic  seismograms  to  the  data.  The 
frequently  stated  justification  for  this  is  that  the  latter  parts  of  the  body  waves  are  considered  to  be 
scattered  waves  with  little  relationship  to  the  source  processes.  This  explanation  is  quite  plausible, 
but  the  length  of  such  windows  in  practice  is  still  determined  quite  arbitrarily.  Some  studies  that 
use  only  about  a  dozen  long  period  waveforms  also  claim  to  have  determined  source  depth  to  an 
accuracy  of  a  few  kilometers  and  having  seen  indications  of  source  finiteness.  Since  the 
wavelengths  of  the  seismic  waves  at  the  frequencies  where  the  waveforms  seem  to  fit  are  much 
larger  than  the  depth  or  the  physical  dimensions  of  the  source  and  because  the  fits  are  not  perfect  at 
any  frequency  such  claims  seem  to  be  poorly  supported.  Apparent  internal  contradictions  and 
inconsistencies  of  this  kind  indicate  a  need  for  some  quantitative  justification  of  the  results  and  the 
procedures. 

Although  many  of  the  apparent  problems  associated  with  time  domain  waveform  inversion 
techniques  using  only  long  period  data  may  appear,  at  a  first  glance,  to  have  little  relevance  to 
nuclear  monitoring,  the  same  kinds  of  problems  need  to  be  solved  in  nuclear  seismology  as  well. 
The  types  of  generic  wavelength-resolution  problems  we  shall  discuss  are  prominent  in  resolving 
source  mechanisms  of  compound  events  using  surface  waves  and  long  period  body  waves,  thus, 
affecting  the  yield  estimation  from  Ms  (Patton,  1988;  Cohee  and  Lay,  1988).  The  techniques 
widely  used  for  source  inversion  from  short  period  data  more  relevant  to  nuclear  monitoring  are  the 
same,  only  the  frequency  ranges  are  different.  For  instance,  particulars  of  the  sources  such  as  pP, 
spall  and  strain  release  that  could  affect  yield  estimation  could  be  derived  better  by  appropriate 
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analyses  of  the  spatio-temporal  behavior  of  explosion  and  small  earthquake  sources.  Admittedly, 
for  such  studies  short  period  or  broadband  data  will  have  to  be  emphasized  more  and  the  analysis 
of  such  data  will  be  even  more  difficult  because  of  the  site  effects.  Our  ultimate  goal  is  to  apply  the 
methods  we  develop  to  explosions  and  small  earthquakes  on  a  large  scale,  but  we  need  to  analyze 
the  less  complex  problem  of  inverting  long  period  data  first.  Another  motivation  for  this  work  is 
that  while  time  domain  modeling  and  waveform  analysis  as  practiced  today  is  a  laborious, 
interactive  and  subjective  trial-and-error  process,  our  procedures  can  be  automated  relatively  easily 
and  may  form  parts  of  future  systems  for  performing  routine  discrimination  and  yield  estimation 
tasks. 


It  seems  certain  that  the  conclusions  of  many  time-domain  studies  could  stand  up  to 
rigorous  analyses  of  errors  and  resolution.  Nevertheless,  many  probably  would  not  and  the 
widespread  use  of  such  techniques  makes  it  desirable  to  develop  more  quantitative  procedures  for 
assessing  the  reliability  and  stability  of  solutions.  If  the  waveform  fits  are  perfect  (and  the  model 
to  be  fitted  is  known)  it  makes  no  difference  whether  the  fitting  is  done  in  the  frequency  or  time 
domain,  but  since  it  is  not  possible  to  fit  any  observed  data  perfecdy,  it  matters  how  well  the  data 
matches  the  predictions  from  the  models  in  either  domain.  This  paper  will  address  some  aspects  of 
this  problem.  Formulated  in  the  frequency  domain  the  inversion  problem  is  completely  frequency- 
separable  and  linear  as  seen  by  rewriting  the  well  known  equations  from  Aki  and  Richards  (1980, 
p.  5r3) 


(  X,  tu )  =  (to ) .  {0})  +  N^{0})  (12) 

where  u  are  displacements,  M  are  components  of  moment  tensors,  G  are  derivatives  of  the  Green's 
functions  and  N  are  noises.  The  problem  in  the  frequency  domain  also  becomes  much  simpler  and 
less  expensive  to  solve  computationally  (Olson  r.nd  Anderson,  1988).  In  frequency  domain 
inversion,  instead  of  taking  the  time  domain  trace  amplitudes  as  data  we  could  fit  the  amplitudes 
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and  phases  of  the  various  frequency  components  at  the  various  recording  stations.  It  is  also  much 
simpler  to  characterize  site  distortion,  S/N  ratios  and  attenuation  in  the  frequency  domain  and 
optimize  an  inversion  process  accordingly.  It  also  follows,  by  simply  inspecting  equation  (12)  that 
information  on  fine  details  found  in  the  space-time  structure  of  seismic  sources  cannot  be  recovered 
from  the  data  unless  the  quality  of  the  synthetic  data  fits  is  acceptable  at  frequencies  which 
determine  the  essential  features  of  the  solution  in  time  or  spatial  wavelength.  Thus,  fitting  the  data 
in  some  frequency  range  does  not  define  the  source  spectrum  at  other  frequencies  although  non¬ 
linear  constraints  often  based  on  plausible  assumptions  or  a  priori  knowledge  about  the  source  may 
influence  the  solution  obtained  at  frequencies  not  fitted.  Using  the  linear  relationship  in  the 
equation  above,  it  is  clear  that  we  can  backpropagate,  in  the  frequency  domain,  any  fitting  errors  to 
assess  the  stability  of  a  source  inversion  solution,  regardless  of  how  complex  the  original  inversion 
process  has  been.  Such  backpropagation  can  provide  the  a  posteriori  error  and  resolution  limits  of 
the  time  domain  source  parameter  estimation  procedures. 

With  regards  to  the  time  domain  inversion  the  following  questions  arise:  Does  the 
existence  of  "good"  visual  fits  from  time  domain  waveform  analyses  as  presently  practiced  ensure 
the  validity  of  deductions  with  respect  to  source  depth,  time  functions  and  spatial  details  of  sources 
as  frequently  presented?  How  well  do  we  have  to  fit  seismic  waveforms? 

In  the  following,  we  shall  try  to  answer  some  questions  of  this  kind  by  examining  how 
well  the  waveform  fitting  is  done  in  common  practice  as  a  function  of  frequency  and  introducing 
some  statistical  tests  to  test  the  stability  of  inversions  and  examine  the  spatial  resolution.  In  this 
paper,  we  shall  focus  on  the  relatively  simple  problem  of  inverting  the  waveforms  of  long  period 
body  waves.  The  main  reason  for  this  is  that  array  studies  have  shown  that  short  period 
waveforms  vary  so  much  over  short  distances  that  fitting  their  waveforms  in  detail  is  hardly 
worthwhile  unless  we  estimate  waveform  distortion  due  to  near-receiver  geology  and  make 
corrections  for  them.  This  is  a  complex  procedure  that  may  make  use  of  signal  frequency 
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components  up  to  4-5  Hz  (Filson  and  Frasier,  1972;  Der  et  al,  1987;  Shumway  and  Der,  1985). 
Path  and  site  related  distortions  of  waveforms  can  place  an  upper  frequency  limit  on  the 
frequencies  to  be  utilized  directly,  i.e.,  without  site  correction,  in  waveform  inversion.  This  limit 
seems  to  be  below  1  Hz,  but  probably  it  is  not  below  0.3  Hz,  where  the  inter- sensor  coherences 
seem  to  be  high  over  large  arrays  such  as  NORSAR.  At  large  arrays  which  were  placed  in 
tectonically  stable  areas  (EKA,  YKA)  the  waveform  similarity  is  still  excellent  at  1  Hz.  The 
statistical  analysis  of  waveform  variability  is  clearly  beyond  the  scope  of  this  paper  and  should  be 
the  subject  of  future  research. 

The  prevalence  of  site  distortion  and  near-source  scattering  is  often  given  as  the  reason  for 
de-emphasizing  the  high  frequency  parts  of  seismic  signals,  often  above  0.15  Hz,  in  long  period 
waveform  inversion  studies  as  we  shall  see  below.  Nevenheless,  combinations  of  long  and  short 
period  waveforms  are  often  fitted  in  time  domain  studies.  This  is  indicative  of  the  fact  that  no 
general  agreement  exists  as  to  the  upper  frequency  limit  above  which  the  data  cannot  be  used  in 
source  studies.  Given  the  observations  quoted  above,  this  requires  further  justification  based  on 
the  study  of  site-related  waveform  variability.  The  reasons  for  choosing  a  frequency  range  for  data 
inversions  do  not  affect  our  conclusions  below;  the  consequences  of  such  practices  in  decreasing 
the  resolution  of  the  source  inversion  results  will  be  the  same.  The  need  for  some  quantification  of 
the  errors  in  such  procedures  has  motivated  some  researchers  to  develop  formal  procedures  for 
error  estimation  (Tichelaar  and  Ruff,  1989).  These  developments  are  important  though  they  are 
still  restricted  to  time  domain  inversion. 

The  purpose  of  this  paper  is  not  to  present  detailed  criticisms  of  individual  studies  in  the 
literature,  but  after  describing  some  generic  problems  with  the  present  approaches  to  the  problem, 
to  outline  frequency  domain  procedures  for  estimating  the  spatial  resolution  of  typical  long  period 
body  wave  inversions  and  show  how  to  test  the  stability  of  solutions  and  estimate  confidence  limits 
for  the  estimates  of  parameters  with  a  given  quality  of  fit  specified  as  a  function  of  frequency. 
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TABLE  V:  Sources  of  Waveform  Synthetic  Pairs  for  Coherence  Analyses 


1)  Goff,  J.A.,  Bergman,  E.A.  and  S.C.  Solomon  (1987).  Earthquake  source  mechanisms  and 
transform  fault  tectonics  in  the  Gulf  of  California.  J.  Geophys.  Res.,  92,  10,485-10,510. 

2)  House,  L.S.  and  K.H.  Jacob  (1983).  Earthquakes,  plate  subduction,  and  stress  reversals  in 
the  eastern  Aleutian  Arc.  J,  Geophys.  Res.,  88,  9347-9373. 

3)  Huang,  P.Y.,  Solomon,  S.C.,  Bergman,  E.A.  and  J.L.  Nabelek  (1986).  Focal  depths  of  mid- 
Atlantic  Ridge  earthquakes  from  body  waveform  inversion.  J.  Geophys.  Res.,  91,579-598. 

4)  Huang,  P.Y.  and  S.C.  Solomon  (1987).  Centroid  depth  and  mechanism  of  mid-ocean  ridge 
earthquakes  in  the  Indian  Ocean,  Gulf  of  Aden,  and  Red  Sea.  J.  Geophys.  Res.,  92,  1361-1382. 

5)  Molnar,  P.  and  W-P.  Chen  (191983).  Focal  depths  and  fault  plane  solutions  of  earthquakes 
under  the  Tibetan  Plateau.  J.  Geophys.  Res.,  88,  1180-1196. 

6)  Nabelek,  J.  (1985).  Geometry  and  mechanism  of  the  1980  El  Asnam,  Algeria,  earthquake 
from  inversion  of  teleseismic  body  waves  and  comparison  with  field  observations.  J.  Geophys. 
Res.,  90,  12713-12728. 

7)  Pelayo,  A.M.  and  D.A.  Wiens  (1989).  Seismotectonics  and  relative  plate  motions  in  the  Scotia 
Sea  region.  J.  Geophys.  Res.,  94,  7293-7320. 

8)  Peppopane,  S.K.  and  S.G.  Wesnousky  (1989).  Large  earthquakes  and  crustal  deformation 
near  Taiwan.  J.  Geophys.  Res.,  94,  7250-7264. 

Disclaimer:  Selection  of  these  papers  for  assessing  the  average  frequency 

characteristics  of  waveform  fits  does  not  necessarily  imply  criticism  of  the 
conclusions  of  these  papers. 


Quantification  of  Time  Domain  Waveform  Fi»s 

To  judge  the  effectiveness  of  waveform  fitting  methods,  we  attempted  to  assess,  in  a 
quantitative  way,  the  quality  of  wavefoim  fits  generally  considered  acceptable.  In  order  to  do  this, 
we  have  taken  21  synthetic  and  data  waveform  pairs  randomly  '^elected  from  the  literature  (a  list  of 
the.se  publications  is  given  in  T^ble  V),  enlarged  and  digitized  them  on  a  digitizing  table.  These 
samples  were  chosen  to  assess  the  average  quality  of  time  domain  waveform  fits  in  the  literature,  to 
be  used  as  a  point  of  reierence.  We  do  not  irnply  that  any  of  the  critical  remarks  above  necessarily 
apply  to  the  studies  these  waveforms  were  taken  from.  The  waveform  pairs  chosen  are  shown  in 
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Figure  27.  Above  each  pair  the  source  of  each  waveform  pair  is  identified  by  the  reference  number 
in  Table  V.  The  time  domain  correlation  coefficient  follows  a  slash  after  this  identifying  number. 
All  the  waveforms  were  rescaled  to  the  same  time  scab  and  the  positive  maximum  of  the  data 
waveform  was  normalized  to  unity,  while  the  synthetic  retained  the  same  relative  amplitude  with 
respect  to  the  data  in  each  pair.  The  differences  between  the  data  and  the  synthetic  traces  are  clearly 
visible  in  the  original  publications,  and  are  not  products  of  the  digitizing,  i.e.,  these  differences  are 
much  larger  than  any  plausible  digitizing  errors.  Since  all  waveforms  were  plotted  with  an  aspect 
ratio  similar  to  the  original  papers,  they  should  be  easily  recognizable  on  inspection.  The  reader 
can  verify,  by  going  back  to  the  original  studies,  that  we  avoided  the  use  of  the  worst  waveform 
fits,  and  selected  the  better  ones. 

We  computed  the  time  domain  correlation  coefficients  between  data  and  synthetic 
waveforms  in  our  sample 


This  measure,  compounded  over  all  the  waveforms,  is  also  the  most  common  choice  for 
the  judging  the  quality  of  fit  in  time  domain  waveform  inversion  studies  if  any  quantitative 
measures  of  waveform  similarity  are  applied  at  all  (see  for  example  Wallace  et  al,  1981).  This 
parameter  does  not  provide  any  information  about  how  good  the  fits  are  between  the  data  and  the 
synthetic  seismograms  as  functions  of  frequency;  it  merely  provides  a  single  number  characterizing 
the  overall  RMS  fit,  and,  therefore,  is  a  poo.'  measure  for  estimating  spatial  resolution.  Generally, 
it  weights  the  various  frequency  components  proportionally  to  the  cross-power  spectra  of  the  two 
traces,  which  is  an  accidental  product  of  the  instru.nent  respon.se,  source  spectrum  and  the  Q 
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operator,  but  does  not  necessarily  reflect  the  information  content  of  the  signals  vs.  frequency. 
Thus  maximizing  the  time  domain  correlation  coefficients  will  generally  result  in  a  suboptimum 
scheme  for  extracting  information  from  the  signals. 

The  time  domain  correlation  coefficients  for  our  data  sample  average  near  0.9  for  our 
selected  sample;  one  sample  is  over  0.95,  and  only  one  is  below  0.8.  The  visual  similarity  is  also 
"good";  i.e.,  the  chosen  waveforms  fit  as  well  or  better  than  the  average  data- synthetic  pairs  in  the 
literature.  The  synthetic  seismograms  contain  somewhat  less  high  frequency  energy  than  the  data 
in  accordance  with  visual  observations.  Moreover,  the  spectra  are  sharply  peaked  (Figure  28) 
giving  the  data  a  narrowband  character. 

Let  us  now  evaluate  the  data  vs.  synthetic  fits  in  the  frequency  domain  utilizing  an 
alternative  similarity  measure  for  the  same  data  set.  This  measure  is  an  ensemble-averaged 
coherence 
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FIGURE  27:  Waveforms  used  for  computing  the  data-synthetic  coherences  in 
Figure  28.  The  different  traces  are  shown  below  each  pair.  The  two  numbers 
next  to  each  pair  are  the  reference  number  in  Table  II  and  the  time  domain 
correlation  coefficient. 
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(14) 


C(cy)  = 


Ci-(Q))S*(C0) 


Vs  d-  (CD  )  d* (O)  )  ^  S-  (Q)  )  S*  (CO  ) 


that  sums  over  waveform  pairs  i  and  applies  spectral  smoothing  as  denoted  by  the  overbar. 
Plotting  this  as  a  function  of  frequency  in  Figure  29,  it  can  be  seen  that  while  the  fits  near  the  peaks 
of  the  spectra  are  good  (but  by  no  means  perfect)  with  a  coherence  a  little  over  0.9,  this  fit 
decreases  in  quality  rapidly  with  increasing  frequency.  Wave  components  above  five  second 
periods  are  hardly  fitted  at  all.  For  the  coherence  calculations  the  sample  pairs  were  tapered  at  the 
trailing  end  by  a  cosine  taper  on  the  last  1/5  of  the  waveform,  zero  filled  to  the  next  common  power 
of  two  in  the  number  of  samples  prior  to  Fourier  transformation.  The  tapering  actually  does 
increase  the  similarity  of  the  waveforms  and  thus  strengthens  our  conclusions  below.  Background 
noise  levels  were  generally  small  compared  to  the  amplitudes  of  the  P  waves  digitized.  By 
equalizing  the  amplitudes  of  the  waveform  pairs  by  the  normalization  described  above,  we  ensured 
that  all  the  sample  pairs  were  considered  roughly  equally  in  the  coherence  analyses  below.  The 
approximate  effective  degrees  of  freedom  for  the  coherence  estimates  shown  were  50  after  spectral 
smoothing. 

Although  we  are  confident  that  our  digitization  is  good  at  least  to  0.3  Hz,  the  coherence  at 
this  frequency  and  below  it  is  essentially  indistinguishable  from  zero  for  this  data  set.  We  shall 
designate  this  coherence  plot  as  representing  the  quality  of  a  "standard  fit"  in  time  domain  modeling 
of  long  period  data.  Clearly,  this  standaid  will  have  a  profound  effect  on  the  resolution  of  the 
time-space  details  of  sources  since  the  best  fits  are  associated  with  frequencies  below  and  near  0.1 
Hz  for  which  the  apparent  wavelengths  are  on  the  order  of  60  km  vertically  and  1 50-250  km 
horizontally,  much  larger  than  the  dimensions  of  the  sources  of  mid-to-large  eanhquakes.  By 
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FIGURE  28:  Averaged  spectra  of  a  sample  of  long  period  data  and  synthetic 
seismogram  pairs  randomly  selected  from  the  literature. 
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FIGURE  29:  Ensemble-averaged  coherences  of  our  sample  of  data-synthetic  pairs 
defining  the  “average  fit”.  The  95%  limit  shows  the  level  below  which  the 
coherences  are  not  significantly  different  from  zero. 
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analyzing  such  waves,  we  are  attempting  to  measure  changes  in  source  depths  on  the  order  of  a 
few  km  and  source  finiteness  on  the  order  of  a  few  tens  of  kilometers. 

We  have  further  verified  our  assessment  of  waveform  fits  by  analyzing  the  differences  of 
the  data  and  synthetic  traces.  Plotting  these,  we  found  that  the  amplitudes  of  the  difference  traces 
suddenly  increased  at  the  arrival  time  of  the  signal  and  remained  at  the  same  level,  about  1/3  or  1/4 
the  signal  amplitude,  throughout  the  waveform  sample  (Figure  1).  Thus,  the  initial  part  of  the 
waveform,  contrary  to  a  false  impression  that  one  might  easily  get  from  visual  inspection  of  the 
waveform  pairs,  is  not  fitted  any  better  than  the  later  part.  The  increase  in  the  amplitude  of  the 
difference  traces  at  the  arrival  time  of  the  signals  shows  that  the  misfit  is  not  due  to  background 
noise.  Moreover,  when  we  compared  the  averaged  spectra  of  the  difference  traces  to  those  of  the 
data  and  synthetic  and  data  traces  themselves  the  spectra  of  the  difference  traces  were  significantly 
below  those  of  the  latter  two  only  in  the  .04-.  14  Hz  range.  This  confirms  our  assessment  from 
coherence  analyses  that  the  average  long  period  data  body  wave  inversions  in  the  literature 
generally  fit  only  a  small  range  of  frequencies  in  a  band  not  exceeding  0.1  Hz.  Thus,  the  fact  that 
the  quality  of  the  waveform  fits  vary  strongly  with  frequency  within  the  pass-band  of  the  long 
period  instrumentation,  which  is  also  quite  evident  from  simple  visual  inspection,  probably  should 
not  require  further  elaboration. 

Statistical  Analysis  of  the  Source  Inversion  Problem 

Let  us  look  now  at  the  space  resolution  obtainable  from  our  typical  long  period  waveform 
fit.  First  of  all,  we  need  some  kind  of  statistical  method  to  find  the  location  of  a  point  source  with 
known  mechanism  in  space  (depth  and  location  within  a  limited  region).  Usually,  there  exists 
some  a  priori  knowledge  of  location  and  depth  from  short  period  data,  the  source  mechanism  may 
be  known  approximately  from  conventional  analyses  of  first  motions.  Out  the  exact  depth  and 
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details  of  source  mechanism  may  be  unknown.  These  could  be  determined  by  finding  the 
maximum  of  the  frequency  domain  statistic 


ISfJ 

F=(N-l) - ^ -  (15) 

// 

which  becomes  a  source  "image"  when  mapped  in  space.  In  this  formula,  N  is  the  number  of 
sensors,  S  is  the  spectral  matrix  of  the  data,  computed  as  d*d,  the  outer  product  of  the  Fourier 
transforms  of  the  data  at  various  sensors,  and/is  the  vector  of  Fourier  components  at  the  various 
sensors  of  the  elementary  source  (synthetic  seismograms)  we  try  to  map  over  space.  Note  that  the 
customary  smoothing  in  constructing  the  spectral  matrix  is  not  an  essential  part  of  this  process  and 
it  can  be,  and  will  be,  omitted  in  the  examples  that  follow.  Thus,  the  method  can  be  applied  to 
short  transients,  such  as  teleseismic  body  waves,  just  as  well.  The  technique  described  above  is 
essentially  an  adaptation  of  a  high  resolution  F-K  algorithm  (Shumway,  1988)  for  finding  the 
position  of  the  elementary  sources  from  waveforms  associated  with  different  slownesses,  instead 
of  the  common  application  of  finding  the  slowness  vectors  of  signals  associated  with  recordings  at 
an  array  with  known  sensor  coordinates.  It  can  be  shown  easily  that  mathematically  the  two 
problems  are  the  same,  since  the  position  and  wavenumber  vectors  are  interchangeable  in  the 
phasor  factors  exp(iitx).  The  additional  complication  caused  by  the  fact  that  the  amplitudes  of  the 
various  frequency  components  are  different  is  taken  into  account  by  normalizing  the  results  of 
beaming  by  dividing  with  ff*.  Thus,  the  problem  of  finding  the  source  position  in  a  limited 
source  region  is  analogous  to  finding  the  wavenumber  vectors  of  noisy  signals  observed  at  an 
array  of  limited  spatial  extent.  Since  in  our  application  the  wavelengths  of  the  signals  are  generally 
larger  than  the  size  of  the  source  region  and  there  is  always  some  mismatch  noise  the  problem  is  an 
extremely  difficult  one.  Therefore,  we  must  expect  serious  problems  due  to  trade-offs,  lack  of 
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resolution  and  sidelobes  in  the  process.  These  problems  are  present  in  the  time  domain  approach 
to  the  problem  as  well,  although  this  may  not  be  obvious. 

Similar  ideas  pertaining  to  finding  the  positions  of  signal  sources  in  complex  waveguides 
were  advanced  by  Baggeroer  et  al  (1988),  Byrne  et  al  (1990)  and  Steele  and  Byme  (1990),  among 
many  others,  in  underwater  acoustics  for  locating  submarines,  and  Goldstein  and  Archuleta  (1987) 
in  seismology  using  the  Capon  and  MUSIC  (Capon,  1969;  Schmidt,  1981)  algorithms  respectively 
for  studying  seismic  sources.  The  starting  point  for  all  the  algorithms  is  the  same,  we  have 
waveform  data  and  some  algorithms  for  computing  Green’s  functions.  The  goal  is  to  estimate  the 
spatial  distribution  of  sources  and  their  time  functions. 

Our  statistic  is  simply  one  of  the  many  possible  algorithms  for  examining  the  spatial 
resolution:  we  chose  it  because  it  is  intuitively  and  computationally  simple,  does  not  require  matrix 
inversion  or  eigenanalysis,  and  has  a  tractable  probability  distribution.  The  F  statistic  defined  by 
equation  (15)  is  essentially  a  ratio  of  the  generalized  beam  power,  the  expression  in  the  numerator, 
and  the  residual  power  or  mismatch  noise  in  the  denominator  (a  S/N  ratio,  not  to  be  confused  with 
the  conventional  definition  where  the  noise  is  commonly  defined  as  ambient  noise).  This  statistic 
for  a  single  frequency  has  an  approximate  non-central  F  distribution  at  individual  frequencies  with 
2  and  2(N-1)  degrees  of  freedom  for  N  input  waveforms  and  the  approximate  associated  non¬ 
centrality  parameter  is  twice  the  estimated  S/N  power  ratio  i.e.,  2FI(N-1)  (Shumway,  1988;  Hogg 
and  Craig,  1978). 

Particulars  of  this  statistic,  its  relationship  to  Capon’s  (1969)  estimator  are  given  in 
Appendix  A  of  this  report.  A  simpler  alternative  algorithm  based  on  a  c^  statistic,  yet  untested  is 
described  in  Appendix  B.  If  we  stack  the  F  statistics  obtained  at  various  frequencies  over  the 
source  coordinates  the  probability  density  function  of  the  result  is  approximately  normal  with  mean 
Zrrii  and  variance  with  E  Ci^the  sums  are  being  over  all  the  frequencies  used,  and  m  i  and  Oi  are 
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the  mean  and  the  standard  deviation  estimates  of  the  individual  non-central  F  distributions  for  a 
single  frequency.  The  degree  of  misfit  will  affect  the  mean  and  the  standard  deviation  of  the  F 
statistic  at  each  frequency,  thus  providing  a  way  to  obtain  a  posteriori  error  estimates  automatically. 
The  method  above  can  also  be  considered  a  backpropagation  algorithm  in  which  the  phase  shifts 
due  to  propagation  are  compensated  for  and  we  are  testing  for  the  phase  alignment  of  the  various 
Fourier  components  at  various  source  coordinates.  Note  that  in  the  algorithm  any  common  spectral 
factors,  such  as  those  associated  with  instrument  responses  and  source  spectra,  divide  out.  Thus, 
it  is  well  suited  to  avoid  trade-offs  between  such  factors  and  the  source  depth  and  location  and  to 
investigate  the  latter  separately.  The  mathematical  derivation  for  an  alternative,  yet  untested 
algorithm  based  on  central  statistic  is  given  in  Appendix  B  of  this  report. 

For  testing  the  method  on  a  simple  case,  we  shall  assume  the  source  crust  to  be  a  half-space 
so  that  the  generalized  beamsteer  vector  becomes 

(  ^  1 

(<y )  =  (16) 

where  the  subscripts  j=l,2,3  refer  to  the  the  phases  P,  pP  and  sP  respectively  and  A:,- is  a  unit 
vector  appropriate  to  a  station  /,  jc  is  a  position  vector  in  the  source  region,  the  Cj  are  phase 
velocities  at  a  given  station,  and  the  tij  are  the  travel  times  of  the  P,  pP  and  sP  phases.  The  source 
"image"  is  a  mapping  of  F  over  the  sourxo  region  in  depth  and  horizontal  coordinates.  Depth 
dependence  enters  into  this  formulation  through  the  delay  times  of  the  depth  phases  and  the 
exponential  factors  are  determined  by  the  horizontal  coordinates  of  the  source  relative  to  some 
reference  point  and  by  the  phase  slownesses.  We  assume  for  the  sake  of  convenience,  that  the 
Green's  functions  are  perfectly  known,  which  is  generally  not  the  case.  The  common  problem  of 
poorly  defined  Green's  functions,  which  makes  the  inverse  problem  even  less  tractable,  is  not 
addressed  in  this  paper. 
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Having  specified  a  statistically  tractable  quantitative  way  to  image  a  source,  we  applied  this 
algorithm  to  a  typical  sample  problem,  that  of  decomposing  a  compound  source  with  known 
mechanism  using  P  wave  data  from  10  stations,  azimuthally  well  distributed  around  this  source 
with  the  source-to-station  azimuths  and  slownesses  listed  in  Table  VI.  The  P  waves  contain  only 
the  P,  pP  and  sP  arrivals.  Our  compound  source  source  in  this  test  consisted  of  two  double 
couples  of  equal  power  with  identical  mechanisms  (stTike=245°,  dip=30°  and  rake  =120°)  the  first 
situated  at  a  depth  of  3  km  near  the  center  of  the  search  grid,  and  the  other  displaced  to  the 
northeast  by  28  km  at  the  depth  of  9  km  in  an  elastic  half-space  with  6.5  km/sec  compressional  and 
3.8  km/sec  shear  velocity.  This  is  about  the  simplest  compound  source  inversion  problem  there  is, 
and  it  is  representative  of  the  scale  of  a  medium-sized  earthquake. 

In  Figure  30,  we  show  images  of  this  source  (or  in  the  parlance  of  underwater  acoustics 
ambiguity  surfaces),  i.e.,  plots  of  F  located  within  a  coarse  search  grid  with  5  km  increments,  at 
different  depths.  Normally,  of  course,  one  would  not  search  an  area  of  250  X  250  km,  such  as  the 
one  shown,  but  would  look  for  maxima  of  some  statistic  in  a  smaller  source  area  or  a  fault  plane. 
Nevertheless,  presenting  such  images  should  give  a  good  idea  of  the  spatial  resolution  provided  by 
the  long  period  waves  we  utilize  here.  For  the  result  shown  in  Figure  30,  we  have  added 
gaussian-distributed  random  complex  numbers  to  the  original /  vectors  (for  a  source  at  6  km  depth) 
to  correspond  to  a  0.9  coherence  (correlation  coefficient)  value  between  the  simulated  "data"  and 
the  synthetics  (the/’s)  and  used  the  frequencies  between  0.02  and  0.15  Hz,  stepped  by  0.0125 
Hz,  with  equal  weight  to  define  the  frequency  range  utilized  by  the  "average  fit."  The  statistic 
gives  a  maximum  that  is  located  between  the  two  sources  and  there  is  very  little  indication,  if  any, 
of  the  presence  of  two  sources  if  we  are  satisfied  with  this  quality  of  fit.  With  the  low  frequency 
end  of  the  spectrum  fitted  preferentially,  there  is  thus,  a  tendency  for  the  image  to  merge  the  two 
sources  together  at  the  "centroid." 
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TABLE  Vf;  Station  Parameters  for  Numerical  Simulations 


Station  # 

Azimuth  /degrees 

Slowness  sec! km 

1 

0. 

0.08 

2 

0. 

0.05 

3 

60. 

0.06 

4 

120. 

0.06 

5 

120. 

0.05 

6 

180. 

0.08 

7 

180. 

0.05 

8 

240. 

0.06 

9 

300. 

0.05 

10. 

300. 

0.08 
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FIGURE  30:  Perspective  views  of  the  variations  of  the  frequency-compounded  F 
statistics  over  the  frequency  range  of  0.02-0.15  Hz  for  a  pair  of  double-couple 
point  sources  with  equal  power.  Uniform  random  noise  was  added  to  simulate  a 
fit  with  an  average  correlation  coefficient  of  0.9.  The  peak  value  occurs  between 
the  two  sources  and  the  two  sources  were  not  resolved.  The  F  statistics  were 
computed  for  various  assumed  source  depths  and  the  sides  of  the  source  region 
views  are  250  km  long.  The  maximum  value  of  the  statistic  and  its  location  in 
gridpoint  coordinates  from  the  left  lower  corner  (5  km  spacing  was  used  in  both  x 
and  y)  are  given  beside  each  view.  The  gridpoint  x-y  coordinates  of  the  two 
synthetic  sources  were  at  27,  27  and  31,  31  and  the  depths  were  at  3  and  9  km 
respectively. 
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FIGURE  31:  Perspective  views  of  the  variations  of  the  frequency-compounded  F 
statistics  over  the  frequency  range  of  0.02-0.35  Hz  for  a  pair  of  double-couple 
point  sources  the  F  values  were  computed  for  various  assumed  source  depths. 
Random  noise  with  16%  of  the  signal  was  added  throughout  the  band.  The  peak 
values  occur  nearer  to  the  two  sources  and  the  two  sources  were  marginally 
resolved  (for  additional  details  of  this  display  see  the  caption  of  Figure  4). 
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We  next  investigated  how  image  could  be  improved  if  we  extended  the  frequency  range  of 
fitting  to  0.02-0.35  Hz.  We  added  16%  noise  throughout  this  band,  thus  doing  away  with  the 
"standard  fit"  and  attempting  to  fit  the  high  frequencies  as  well  as  the  low  frequencies.  As  Figure 
31  indicates,  we  now  see  two  maxima  at  two  different  depths  corresponding  to  the  two  sources. 
The  locations  are  still  closer  to  each  other  than  the  actual  ones,  and  we  still  have  a  sizable  peak 
between  the  two  source  locations  at  depth  where  we  put  no  sources. 

Finally,  we  decided  to  discard  the  lower  range  of  frequencies  altogether  and  did  the  fitting 
in  the  0.15-0.35  Hz  range  with  16%  noise  added.  We  see  in  Figure  32  that  the  two  sources  well 
resolved  with  the  values  of  the  statistic  between  the  two  sources  much  diminished,  but  the 
tendency  of  the  two  sources  to  merge  is  still  evident  by  the  biases  in  the  positions,  the  peaks  being 
closer  together  than  the  originally  specified  source  locations. 

These  results  spell  trouble  for  the  studies  of  moderate-to-large  earthquakes,  based  on  long 
period  data  only,  since  it  appears  that  fitting  the  data  does  not  allow  us  to  resolve  details  of  fault 
motion.  These  results  do  not  imply  that  one  could  not  outline,  with  the  prevalent  quality  of 
waveform  fitting,  the  faulting  history  of  larger  events  with  source  dimensions  considerably  larger 
than  the  location  and  depth  uncertainties  found  in  our  simulations.  Nevertheless,  even  in  such 
cases  not  much  detail  in  the  faulting  process  could  be  derived  from  the  data.  The  most  that  can  be 
achieved  using  only  the  low  frequency  components  of  the  data  is  some  definition  of  a  "centroid" 
location  and  an  ill  defined  associated  source  function  which  may  indicate  by  its  long  duration  that 
the  source  had  a  significant  spatial  dimension. 

One  may  ask  whether  our  choice  of  the  test  statistic  could  be  the  cause  of  poor  resolution 
and  the  application  of  the  more  conventional  statistic,  the  time  domain  correlation  coefficient  could 
have  yielded  better  resolution.  In  Figure  33,  we  show  perspective  plots  of  the  expected  correlation 
coefficient  for  various  locations  and  at  a  few  depths  obtained  by  fitting  a  single  double  couple  as  a 
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FIGURE  32:  Perspective  views  of  the  variations  of  the  frequency-compounded  F 
statistics  over  the  frequency  range  of  0.15-0.35  Hz,  thus  deleting  the  low 
frequencies,  for  a  pair  of  double-couple  point  sources  the  F  values  were 
computed  for  various  assumed  source  depths.  Random  noise  with  16%  of  the 
signal  was  added  throughout  the  band.  The  peak  values  occur  near  to  the  two 
sources  and  the  two  sources  were  well  resolved  (for  additional  details  of  this 
display  see  the  caption  of  Figure  4). 
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FIGURE  33:  Perspective  views  of  the  variations  of  the  positive  parts  of  the 
simulated  site-averaged  time  domain  data  vs.  probe  vector  (data-synthetic) 
correlation  coefficient  over  the  frequency  range  of  0.02-0.15  Hz  for  a  pair  of 
double-couple  point  sources  the  F  values  were  computed  for  various  assumed 
source  depths.  The  negative  parts  of  the  correlation  coefficient  were  set  to  zero 
arbitrarily,  since  only  peak  positive  values  are  of  interest.  Random  noise  was 
added  throughout  the  band  to  give  a  peak  value  near  0.9.  The  two  sources  were 
not  resolved  (for  additional  details  of  this  display  see  the  caption  of  Figure  4). 
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source.  The  time  domain  correlation  coefficient  was  approximated  by  the  frequency  domain 
equivalent  formula  (we  have  set  up  the  whole  simulation  algorithm  in  the  frequency  domain)  for 
zero  lag 


where  Pio))  is  the  power  spectrum  of  the  data  (those  of  the  synthetics  are  very  similar)  shown  in 
Figure  2,  the  sums  arc  over  frequency,  the  d  are  the  Fourier  transforms  of  the  synthetic  data  vector 
(sums  of  direct  P  and  surface  reflections  for  all  the  sources),  and  the /are  the  probe  vectors  (in  our 
case,  the  Fourier  transforms  of  the  expression  in  equation  5  for  each  sensor).  In  order  to  resolve 
the  multiplicity  of  this  source,  one  would  expect  that  the  correlation  coefficient  should  be  small 
enough  to  alert  us  that  this  is  not  a  simple  source,  moreover,  we  would  hope  that  the  correlation 
coefficient  would  show  two  distinct  peaks,  each  at  the  proper  source  position.  None  of  these 
expectations  are  met,  however.  The  maximum  values  of  the  correlation  coefficient  are  closer  to  the 
deeper  source,  for  this  realization  there  are  two  identical  maxima  of  0.897  at  6  and  9  km  depth  and 
the  value  at  3  km  depth  is  smaller,  the  locations  do  not  correspond  well  to  the  actual  positions  of 
the  sources.  Accepting  the  standard  of  0.9  correlation  coefficient  one  would  end  the  search  and 
announce  an  acceptable  fit  for  a  simple  double  couple  source  at  these  depths  and  either  of  these 
locations.  The  correlation  coefficient  gives  no  clue  about  the  presence  of  the  shallower  source, 
there  is  no  double  maximum.  Note  the  large  sidelobes  of  the  correlation  coefficient  which  are  due 
to  the  limited  range  of  frequencies  applied  and  their  long  wavelengths.  Incidentally,  both  the  F 
statistic  and  the  correlation  coefficient  are  expected  to  behave  similarly  when  the  frequencies  used 
are  limited  to  a  narrow  range.  The  generalized  beam  on  which  the  F  statistic  depends  is  maximized 
at  the  same  spatial  points  where  the  correlation  coefficient  is  maximized,  and  as  we  have  shown 
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above  none  of  them  resolve  the  two  sources  from  the  information  below  0.15  Hz.  On  the  other 
hand,  when  the  frequency  range  is  wider  the  F  statistic  can  keep  track  of  the  spatial  resolution  by 
estimating  and  compounding  the  same  for  all  frequency  components,  while  the  correlation 
coefficient  does  not  convey  the  same  information  since  the  same  correlation  coefficient  can  be 
obtained  by  mismatching  different  combinations  of  frequency  components. 

Our  results  for  the  simple  problem  of  resolving  two  not  very  closely  spaced  sources  of 
equal  strength  in  a  half-space  imply  that  the  prevalent  fitting  standards  and  frequency  components 
utilized  are  just  on  the  verge  of  indeterminacy  for  resolving  spatial  details  of  seismic  sources  of 
moderate-to-large  events  from  typical  long  period  body  wave  waveform  inversions.  The 
indeterminacy  often  manifests  itself  in  the  fact  that  in  many  studies,  if  the  fits  are  deemed 
unsatisfactory,  adding  a  few  additional  parameters  or  complications  to  the  source  model  will 
generally  produce  good  waveform  fits.  If  such  problems  were  overdetermined,  i.e.,  there  was 
enough  redundancy  in  the  data  set,  then  it  would  be  more  difficult  to  find  solutions  and  the  inability 
to  find  a  solution  might  sometimes  point  to  the  inadequacy  of  the  models. 

The  preference  for  fitting  the  low  frequency  end  of  the  spectrum  is  based  on  the  perception 
that  low  frequency  components  of  the  signal  appear  to  be  more  stable.  One  could  argue,  however, 
that  since  in  studies  of  source  spatial  characteristics  we  essentially  try  to  measure  small  time  delays 
among  subevents  and  various  other  arrivals  (such  as  pP  and  sP)  the  phases  of  low  frequency 
components  (and  the  associated  waveforms)  will  be  less  sensitive  to  variations  in  these  times  than 
high  frequency  components.  Thus,  at  least  sonte  of  this  stability  may  be  illusory. 

Depth  Estimation  Error  Limits  for  the  El  Golfo  Earthquake 

Our  choice  fell  on  this  event  since  the  El  Golfo  earthquake  was  reported  to  be  a  simple 
point  source,  a  strike  slip  event  found  to  be  at  10  km  depth  by  Ebel  et  al  (1978)  and  thus,  similar  to 
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the  synthetic  case  analyzed  above.  The  source  time  function  was  estimated  by  Ebel  et  al,  to  be  a 
simple  triangle  of  4  sec  duration  and  the  t*  along  the  paths  to  the  receiving  stations  recording  on 
WWSSN  15-100  LP  instruments  was  assumed  to  be  1.3  sec.  We  adopted  these  parameters  as 
well  as  the  fault  orientation  and  the  crustal  model  given  in  the  same  paper  and  used  formula  (4)  in 
an  attempt  to  verify  the  depth  estimate  from  the  P  wave  data  by  our  imaging  method. 

The  Green's  functions  (f  vectors)  were  calculated  by  using  an  adaptation  of  the  formulation 
of  body  wave  radiation  from  a  point  source  given  by  Douglas  et  al  (1971)  generalized  to  a  flat¬ 
layered  source  crust  model.  The  synthetic  seismograms  generated  are  very  close  to  the  ones  given 
in  the  paper  by  Ebel  et  al  (1978).  We  attribute  some  small,  occasional  differences  to  ti.e  fact  that 
while  our  program,  based  on  propagator  matrices,  automatically  includes  all  rays  in  the  source 
crust,  the  program  used  by  Ebel  et  al  needs  the  specification  of  all  rays  desired  (Langston,  1976). 
Thus,  some  less  important  rays  may  have  been  omitted  in  the  calculations  of  Ebel  et  al.  In  any 
case,  the  synthetic  waveforms  we  have  generated  are  matching  the  data  just  as  well  as  those  of  Ebel 
et  al  (1978)  and  are  very  similar  to  those  used  in  that  study. 

The  long  period  P  waveforms  to  be  processed  were  obtained  by  digitizing  the  enlarged 
seismograms  from  the  paper  on  a  high  resolution  digitizing  tablet.  We  have  matched  a  37  sec  long 
segment  of  the  waveforms,  similar  to  the  original  study  of  Ebel  et  al  (1978),  and  tapered  the 
trailing  edge  with  a  smooth  cosine  taper.  This  essentially  utilized  the  first  1-1/2  cycle  of  the 
signals,  similar  to  the  range  of  the  visual  fitting  in  the  original  paper.  Compared  to  the  sizes  of  the 
mismatches  between  synthetic  and  data  the  digitizing  errors  and  any  trace  distortions  are  very  small 
and  should  not  affect  the  our  conclusions.  In  order  to  eliminate  the  absolute  times  as  a  factor  in  the 
matching  procedure  (had  we  known  the  absolute  times,  properly  corrected  for  mantle 
heterogeneity,  we  could  have  obtained  the  depth  from  those  alone),  we  lined  up  both  the  time 
domain  representations  of  the/s  and  the  data  on  the  first  P  arrival. 
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We  have  made  two  test  runs  while  varying  the  source  depth  an  fitting  the  spectral  structures 
in  the  .02  to  .15  Hz  band  using  all  frequencies  with  equal  weight.  In  the  first  test,  we  matched  the 
synthetic  seismograms  for  10  km  depth  with  source  crustal  responses  appropriate  to  all  depths  and 
added  5%  random  white  noise  (to  avoid  a  singularity  at  10  km  depth)  and  in  the  second  we 
matched  the  data  with  the  source  crustal  responses  appropriate  for  different  source  depths;  for  the 
results  of  the  latter  we  also  computed  our  admittedly  crude  error  bounds  analytically,  and  assumed 
that  the  compound  F  statistics  have  an  approximate  normal  distribution.  The  means  and  variances 
of  this  normal  distribution  were  computed  by  summing  for  all  the  various  frequency  components 
the  means  and  variances  of  the  individual  non-central  F  distributions,  appropriate  to  the  quality  of 
fit  (S/N  ratios)  implied  by  the  estimated  values  of  F  statistics  (Hogg  and  Craig  1978,  p.  288).  The 
means  and  variances  of  these  individual  non-central  F  distributions  were  computed  by  utilizing  the 
approximate  formulas  given  by  Mudholkar,  Chaubey  and  Lin  (1976).  Thus,  for  estimating  the  the 
error  bounds  the  quality  of  fit  is  automatically  taken  into  account  by  our  algorithm. 

As  one  would  expect,  the  computation  on  the  synthetic  seismogram  for  10  km  depth 
resulted  in  an  unambiguous  maximum  of  the  statistic  at  10  km  depth  with  narrow  95%  error 
bounds  relative  to  the  F  value  at  the  peak  (Figure  34a).  The  data,  on  the  other  hand,  when 
subjected  to  the  same  analysis  procedure  could  be  reconciled  with  practically  any  depth  shallower 
than  16  km  (Figure  34b)  and  had  much  smaller  F  values.  There  are  two  peaks  at  5  and  14  km 
depth,  but  they  are  not  very  significant.  Some  of  the  depths  in  the  acceptable  range  are 
considerably  outside  of  the  2  km  accuracy  limits  claimed  by  Ebel  at  al  (1978)  for  their  solution 
which,  at  10  km,  is  not  inconsistent  with  our  depth  limits.  Theoretically,  if  the  synthetic 
seismograms  and  the  data  are  similar  enough  they  should  give  the  same  solution.  It  is  shown 
above  that  the  data  and  the  synthetic  seismograms  for  the  same  event,  although  their  waveforms  are 
similar,  give  substantially  different  conclusions  with  regards  to  source  depth  which  could  have 
caused  problems  if  we  had  tried  to  analyze  the  fault  motion  in  any  detail  on  the  fault  plane.  We 
have  essentially  shown  in  Figure  34a  that  for  a  similar  problem  the  method  should  have  given  the 
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same  answer  had  the  waveforms  been  more  similar.  Despite  the  fact  that  the  original  El  Golfo 
earthquake  study  also  utilized  some  additional  data,  a  limited  amount  of  SH  and  surface  waves,  we 
do  not  believe  that  the  inclusion  of  these  would  have  changed  our  results  or  conclusions.  We  agree 
with  Ebel  et  al  (1978)  that  the  best  visual  fits  for  P  waves  appear  to  be  near  10  km  in  source  depth. 
We  have  found,  however,  that  this  may  be  deceptive.  In  Figure  35,  we  show  the  data  waveforms 
and  synthetic  seismograms  appropriate  to  various  source  depths  for  the  stations  CAR,  WES,  TRN, 
UME  and  COL.  It  seems  that  the  best  waveform  fit  does  not  occur  at  the  same  depth  for  the 
various  stations.  While  the  best  fit  for  CAR  is  at  10  km,  WES  and  TRN  can  be  reconciled  with  a 
greater  depth  almost  as  well,  while  UME  and  COL  seem  to  fit  the  4  km  synthetic  much  better  than 
than  the  three-peaked  broad  waveforms  for  10  km  depth.  The  net  result  of  such  waveform  fit 
trade-offs  among  stations  gives  an  intuitive  explanation  to  our  results  in  Figure  34.  Although  we 
feel  that  time  domain  waveform  comparisons  are  not  the  best  means  inversion  results  should  be 
judged  by,  we  present  these  waveforms  for  the  benefit  of  those  who  like  to  think  in  the  time 
domain. 

We  used  this  event  to  illustrate  how  confidence  limits  on  source  parameters,  in  our  case  the 
source  depth,  can  be  computed  for  actual  data.  We  were  not  interested  in  the  actual  mechanism  and 
parameterization  of  this  particular  event  which  may  be  refined  by  further  study  that  we  are  not 
planning  to  undertake.  Our  main  point  in  analyzing  this  event  was  also  to  show  that  synthetic  and 
data  waveforms  can  give  different  results  even  if  the  waveforms  are  similar.  Based  on  Monte- 
Carlo  simulations  described  above,  depth  estimation  with  an  accuracy  of  practical  utility  for  events 
of  this  size  and  using  the  frequency  ranges  involved  in  time  domain  studies  of  long  period  body 
waves  would  require  a  fitting  at  least  at  the  .96  coherence  (or  correlation  coefficient)  level.  At  this 
level  of  fitting  accuracy  we  probably  could  not  ignore  the  receiver  crustal  transfer  functions  (as  it  is 
done  presently  in  all  source  inversion  studies).  Thus,  it  appears  that  much  more  work  is  needed 
for  producing  accurate  inversions  than  that  contained  in  most  studies. 
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FIGURE  34:  a)  F  statistics  computed  for  the  set  of  synthetic  seismograms 
appropriate  to  the  El  Golfo  crustal  and  point  source  model  at  10  km  depth  as 
fitted  to  a  wide  range  of  depths.  Five  percent  random  white  noise  was  added  ^ 
avoid  singularity  in  the  solution.  The  dashed  lines  are  approximate  55% 
confidence  limits  to  the  F  values.  The  maximum  unambiguously  defines  a 
of  10  km.  b)  F  statistics  computed  using  the  actual  P  waves  from  the  El  Golfo 
earthquake  and  fitted  to  a  wide  range  of  depths.  The  dashed  lines  are 
approximate  95%  confidence  limits.  The  arrows  indicate  the  large  acceptable 
depth  range  which  is  consistent  with  the  F  values. 
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FIGURE  35;  Selected  data  and  synthetic  seismograms  appropriate  to  various 
source  depths  for  the  El  Golfo  earthquake.  The  “best”  time  domain  waveform  fits 
occur  at  different  depths  at  the  various  stations. 
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Imaging  of  the  1988  Armenian  Earthquake 


This  event  was  described  by  Pacheco  et  al  (1989)  and  was  found  to  be  a  superposition  of 
three  subevents  with  different  source  mechanisms  (Table  2  in  the  paper  by  Pacheco  el  al  1989). 
The  first  of  the  three  subevents  occurred  near  the  center  of  the  causative  fault  break,  it  was 
followed  by  the  second  subevent  located  ESE  from  the  first.  The  third  subevent  occurred  farther 
west  from  the  first  one.  In  this  report,  we  did  not  plan  to  redo  all  the  detailed  analyses  by  Pacheco 
et  al  (1989),  instead  we  only  attempted  to  find  corroborating  evidence  for  their  solution  based  on 
our  methodology.  In  order  to  do  that,  we  have  adopted  the  source  mechanisms  and  depths  for 
their  three  subevents  and  tried  to  verify  the  triple  events  by  matching  each  to  the  P  wave  data  for 
the  eight  broadband  stations  TOL,  COL,  HIA,  BJI,  WMQ,  LZH,  KMI  and  GAR  using  equation 
(15).  The  data  were  obtained  on  a  standard  distribution  tape  for  this  event  from  IRIS  together  with 
the  specification  for  all  the  instrument  responses.  Although  the  three  solutions  are  distinct,  they  do 
not  appear  entirely  dissimilar  at  our  eight  stations. 

Before  applying  our  matching  procedures  we  have,  therefore,  done  some  experiments  on 
the  synthetic  sources  we  generated.  Our  synthetics  were  computed  with  the  source  crust  model  of 
Pacheco  et  al  (1989),  a  t*  of  0.7  sec  that  also  used  by  them,  but  did  not  use  the  trapezoidal  source 
time  functions  applied  in  that  paper  (Pacheco  personal  communication)  since  the  source  time 
function  cancels  in  the  matching  process.  First  of  all,  we  have  matched  each  source  model  to  the 
space-shifted  versions  of  themselves.  This  resulted,  of  course,  in  peaks  placed  at  the  appropriate 
locations  (Figures  36a-c).  Besides  the  peaks  there  are  ridges  running  SW-NE  in  these  images, 
these  artifacts  that  are  related  to  the  fact  that  all  of  our  stations  are  in  the  NW-SE  quadrants  and 
thus,  we  have  poor  control  over  the  inferred  source  locations  in  the  SW-NE  direction.  Matching 
the  three  sources  against  each  other  (Figures  37a-c)  resulted  in  cross-matched  peaks  appropriately 
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FIGURE  37c:  Model  3  against  Model  2.  Nole  that  the  distance  scales  are 
different  on  the  horizontal  (EW)  and  the  vertical  (NS)  axis.  Perspective  views 
are  on  top  and  contour  plots  of  the  same  surface  are  below.  The  sizes  of  grids 
are  127.5  km  x  127.5  km  with  51  subdivisions  in  each  direction. 


placed  on  the  positions  given  by  Pacheco  et  al  (1989).  It  appears  that  sources  1  and  3  match  each 
other  better  than  source  2  against  the  other  two  given  this  staticii  distribution. 

In  the  next  step,  we  have  superposed  the  three  synthetic  sources  and  matched  each  source 
against  this  sum.  We  succeeded  to  recover  the  approximate  locations  of  the  individual  sources, 
although  the  cross-correlation  peaks  also  came  into  play. 

Finally  we  used  the / for  each  of  the  individual  source  type  in  the  matching  process  on  the 
original  data.  Instead  of  deconvolving  instrument  responses  from  each  seismic  trace  of  the  data  or 
constructing  broadband  traces  from  several  instruments  the  components  of  /  vectors  in  equation 
(15)  were  made  to  include  the  appropriate  instrument  responses.  This  approach  should  have  no 
effect  on  the  results  with  regards  to  the  spatial  structure  of  the  sources,  but  has  the  advantage  that 
the  effective  instruments  on  the  various  kinds  of  seismic  recordings  do  not  have  to  be  equalized. 
Matching  the  data  against  the  three  kinds  of  point  sources  (Figures  38  a-c)  results  in  an 
running  ridge  that  faintly  suggests  the  trend  of  the  fault  thought  to  be  responsible  for  the  Armenian 
earthquake  (Pacheco  et  al  1989).  The  peaks  matched  with  subevents  1  and  3  are  located  at  the 
center  and  the  west  end  of  the  fault,  respectively,  as  determined  by  Pacheco  et  al  (1989).  There  is 
no  clear  indication  of  subevent  2.  Since  we  have  used  P  waves  solely  from  only  eight  stations  the 
spatial  resolution  is  poor.  We  plan  to  redo  this  analysis  using  a  much  larger  data  set  approximating 
the  scope  of  the  total  P  and  S  wave  data  used  by  Pacheco  et  al  (1989). 

A  complete  inversion  of  all  the  available  data  from  this  earthquake  is  beyond  the  scope  of 
the  work  on  this  project.  Since  this  is  a  compound  source  a  complete  analysis  would  involve  the 
stripping  out  of  individual  sources  and  test  the  residual  waveforms  for  the  presence  of  other 
sources  taking  into  consideration  of  the  reduced  degrees  of  freedom.  Moreover,  by  estimating  the 
time  functions  of  the  seismic  radiation  from  various  parts  of  the  fault  more  details  about  the  fault 
motion  could  be  revealed. 


115 


Imaging  of  selected  nuclear  explosions  at  Kazakh 


In  order  to  test  our  scheme  of  source  imaging,  we  have  selected  two  Kazakh  events  for 
which  we  had  deconvolutions  at  three  UK  arrays  (Der  et  al,  1987).  Taking  the  deconvolved  traces 
at  these  three  arrays,  we  have  used  the  model  of  an  explosion  source  in  a  half-space  to  design  the / 
vectors  in  the  procedure.  Since  we  have  only  three  waveforms,  complete  imaging  in  the  horizontal 
plane  was  not  practical  since  it  is  not  possible  to  make  a  good  image  with  only  three  azimuths. 
Therefore,  we  have  tried  to  match  the  depth  only.  We  must  point  out  that  while  we  are  trying  to 
match  an  ideal  model  the  actual  deconvolutions  may  not  contain  a  true  pP  in  the  sense  of  an  elastic 
model.  Moreover,  even  if  the  simple  elastic  model  were  applicable,  the  bandwidth  limitations  of 
the  data  may  not  allow  us  to  resolve  the  true  depth.  These  caveats  notwithstanding,  if  we  plot 
Shumway’s  statistics  as  functions  of  source  depth  they  indicate  shallow  source  depths  (Figures  39 
and  40).  Thus,  the  deconvolved  seismograms  best  resemble  an  elastic  compressional  source  at 
shallow  depths;  this  is  more  than  we  can  deduce  from  location  calculations  at  teleseismic  distances. 

It  appears  that  our  imaging  procedure  may  be  a  good  alternative  for  the  discrimination  and 
depth  estimation  schemes  devised  by  Pearce  (1980, 1977)  Pearce  and  Stewart  1989  who  also  tried 
to  match  the  various  depth  phases  with  amplitude  limits  obtained  from  measurements  made  on 
seismic  records.  Our  approach  differs  from  his  in  that  we  can  use  the  whole  seismograms 
automatically  without  manual  measurements,  we  also  utilize  the  polarities  of  phases,  but  without 
trying  to  determine  them  subjectively,  and  that  we  can,  in  addition,  incorporate  any  available 
knowledge  of  crustal  structures,  no  matter  how  complex,  into  the  scheme  by  using  them  to 
compute  the /  vectors.  The  matching  procedure  will  work  in  spite  any  complexities  in  waveforms 
caused  by  non-impulsive  source  time  functions.  Thus,  our  new  standard  method  of  discrimination 
will  require  that  we  line  up  the  P  waves  on  the  first  onset  (the  only  manual  step),  invoke  the  best 
crustal  response  of  the  source  region  and  match  these  against  sets  of/’s  computed  for  explosion 
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FIGURE  39:  Shumway’s  statistic  as  a  function  of  depth  computed  from 
deconvolved  P  wave  seismograms  of  the  December  25,  1975  Kazakh  nuclear 
explosion. 
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FIGURE  40: 
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Shumway’s  statistic  as  a  function  of  depth  computed  from 
wave  seismograms  of  the  June  29,  1977  Kazakh  nuclear 
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sources  at  various  depths  in  that  structure.  In  more  general  crustal  structures,  complexities  not 
explainable  to  simple  depth  phases  will  also  be  present,  and  these  can  also  be  matched.  Good 
matches  (high  F  values)  will  designate  explosions,  poor  matches-various  crustal  earthquakes 
(because  the  mismatches  in  the  amplitudes  an  polarities  of  the  various  arrivals,  mostly  of  the  depth 
phases).  Deconvolution  is  not  needed  for  the  success  of  the  scheme,  because  common  frequency 
factors  cancel  in  the  process;  on  the  other  hand,  corrections  for  unequal  path  attenuation  or 
instrument  responses  will  be  required.  Reduction  of  site  effects  by  array  averaging  is  desirable 
prior  to  this  process,  but  does  not  seem  to  be  absolutely  necessary  if  many  stations  are  being  used. 
Crustal  structures  have  been  the  subject  of  studies  over  many  decades  and  have  been  mapped  all 
over  the  world.  A  database  of  such  stmctures  could  be  a  good  accessory  to  intelligent  monitoring 
systems. 


Analyses  of  surface  waves  from  Kazakh  nuclear  explosions 

We  wish  to  analyze  surface  waves  observed  at  a  set  of  stations  in  a  similar  way,  then  the 
components  of  the  /  vector  become 

fii(0)  =  Pfo) )  exp  \j  Q)k-xl  c- (£0 )  ] 

In  this  case  the  waves  are  dispersive  and  the  phase  velocity  depends  on  frequency,  and 
Pi(0))  is  the  source  radiation  pattern  in  the  direction  of  the  station  i.  The  exponential  again 
describes  horizontal  source  translation  in  laterally  homogeneous  models. 

The  key  problem  in  nuclear  monitoring  is  the  separation  of  isotropic  (explosion)  and 
deviatoric  (double  couple,  strain  release  or  block  motion)  contributions  in  a  .set  of  seismic  surface 
waves.  An  obvious  application  is  yield  estimation  from  surface  waves  when  tectonic  strain  release 
is  present.  We  want  to  reduce  the  effects  of  the  double-couple  component  and  estimate  the  energy 
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in  the  compressional  source  contribution.  For  this  application  the  Capon’s  algorithm  seems  to  be 
the  most  suitable,  since  it  is  essentially  a  formulation  for  the  energy  in  an  output  of  a  constrained 
optimum  filter  designed  to  reduce  the  noise  (double  couple  component)  and  pass  the  signal 
(compressional  component)  undistorted.  The  original  formulation  of  Capon's  F-K  algorithm  has 
the  inverse  roise  spectral  matrix  N  instead  of  that  of  data  S.  In  most  imaging  applications 
(Baggeroer  et  al,  1988)  of  the  total  data  spectral  matrix  containing  contributions  from  both  the 
desired  signal  and  the  undesirable  noise  (strain  release  in  our  case)  is  substituted.  This  may  have 
serious  consequences,  however,  in  the  case  of  mismatch  (Cox,  1973).  In  our  simulations  of 
extracting  the  compressional  contribution  we  shall  use  the  original  formulation  of  Capon  (1969) 
and  will  use  the  spectral  matrix  constructed  only  from  a  priori  estimates  of  the  average  strain 
release  sources  summed  over  a  range  of  mechanisms  {Idi*di)  for  thrusi  events  of  the  type  typical 
for  Kazakh.  According  to  the  literature  of  Capon's  estimator  (Cox,  1973)  the  resulting  algorithm 
that  uses  a  "noise"  spectral  matrix  and  does  not  include  the  desired  compressional  source  will  be 
more  stable  and  much  less  sensitive  to  mismatch.  It  will  be  generally  possible  to  obtain  some 
estimate  of  the  average  strain  release  spectral  matrix  for  a  test  site,  although  this  estimate  need  not 
be  vei7  accurate.  Using  this  in  conjunction  with  an  imperfectly  known  compressional  component, 
we  still  should  be  able  to  reduce  the  contribution  from  the  deviatoric  component  whatever  its  cause 
is. 


Table  VII  shows  the  results  of  such  a  simulation  where  a  constant  compressional  source 
superposed  on  variable  levels  of  a  strain  release  source  (thrust  fault)  was  processed  by  Capon's 
algorithm.  The  "noise"  spectral  matrix  was  constructed  of  a  superposition  of  thrust  fault 
contributions  (that  also  included  the  simulated  fault  used  in  the  input)  the  fault  orientations  of 
which  were  varied  over  a  range  of  five  degrees  to  account  for  the  uncertainties  and  variations  in  the 
strain  release  components.  The  table  shows  that  the  power  in  the  output  (signal  estimate)  did  not 
vary  much,  although  the  ratio  of  the  compressional  and  "strain  release"  components  did.  Only  at 
extremely  large  ratios  can  we  see  a  breakdown  of  the  process.  This  indicates  that  such  a  processor 
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could  work,  despite  the  uncertainties  and  variations  in  the  strain  release  components  as  long  as  we 
have  an  idea  what  we  need  to  eliminate.  The  other  major  issue,  that  of  the  limited  ability  to 
distinguish  various  types  of  shallow  deviatoric  sources  (Patton,  1988)  may  also  be  attacked  by 
similar  signal  processing  methods.  In  the  future,  we  plan  to  investigate  the  associated  resolution 
problem  perhaps  including  long  period  body  waves  as  well  (Cohee  and  Lay,  1988).  In  contrast  to 
the  presently  used  painstaking  procedures  of  fitting  radiation  patterns  to  the  data,  estimating  and 
subtracting  the  strain  release  components,  this  procedure  can  be  automated  and  performed 
instantaneously  after  receiving  the  surface  waves  from  an  event. 


TABLE  VII:  Separation  of  Compressional  and  Deviatoric  Sources  with  Capon's 

Algorithm. 


Relative 

deviatoric 

component 

Output 

amplitude 

0. 

1.0 

0.17 

1.0 

0.333 

0.99 

0.5 

0.99 

1.0 

1.1 

10. 

1.4 

The  only  problem  that  needs  to  be  solved  is  to  obtain  in  practice  good  prior  estimates  of  the 
compressional  and  strain  release  spectral  matrices.  These  can  be  obtained  empirically,  either  by 
finding  events  with  low  strain  release  (small  L^'ve  waves)  or  manipulating  several  events  with 
similar  strain  release  mechanisms,  but  different  degrees  of  strain  release.  Such  procedure  involves 
translating  the  two  events  to  the  same  location  by  phase  shifting  and  differencing  their  spectra  to 
obtain  estimates  of  a  pure  compressional  source.  Such  a  process  would  totally  eliminate  the  need 
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for  complex  "path  corrections.”  It  is  possible,  however,  that  the  method  may  not  work  over  larger 
areas  (exceeding  70  km  in  diameter)  due  to  waveform  distortions  caused  scattering  of  Rayleigh 
waves  (Rivers  and  von  Seggem  1979).  This  would  not  detract  from  the  usefulness  of  such  an 
approach,  since  other  methods  such  as  that  of  Stevens  (1986)  would  be  also  affected  by  the  same 
problem. 

In  view  of  some  data  problems  discovered  and  described  below,  in  this  report  we  set  a  less 
ambitious  goal,  that  of  locating  three  events  relative  to  each  other.  We  assume  that  the  strain 
release  component  is  negligible.  During  this  phase  of  the  work,  we  have  collected  surface  wave 
recordings  of  Kazakh  nuclear  explosions  at  common  stations  using  the  databases  at  the  Center  for 
Seismic  Studies.  We  have  collected  surface  wave  data  for  nine  Kazakh  explosions  which  were 
recorded  at  five  common  stations  ANTO,  GRFO,  MAT,  TATO  and  KONO.  The  dispersive  phase 
shifts  were  estimated  from  the  fundamental  Rayleigh  wave  dispersion  curves  computed  for  the 
Kazakh  crustal  model  given  by  Stevens  (1986).  The  first  test  was  to  verify  that  the  events  gave 
maxima  near  to  their  actual  locations.  This  was  found  to  be  the  case  for  several  of  them.  In  Figure 
41  we  show  the  maps  of  Capon's  statistics  obtained  by  using  the  event  with  the  smallest  tectonic 
release,  lowest  Love  wave  energy,  as  a  master  event  defining  the  /  vector.  This  seems  to  confirm 
the  findings  of  von  Seggem  (1972)  who  was  able  to  locate  events  relative  to  each  other  with 
surface  waves  with  about  a  10  km  accuracy.  However,  some  other  events  did  not  locate  well  at 
all.  Upon  examining  the  data,  we  have  found  that  the  absolute  times  given  in  the  database  were 
grossly  in  error  for  at  least  two  of  the  traces.  Further  inquiries  made  revealed  that  people  who  had 
used  the  data  before  (Richard  Baumstark  for  instance)  also  experienced  problems  with  the  absolute 
times  in  some  of  the  CSS  databases.  While  gross  errors  are  easy  to  notice,  smaller  errors  may  go 
unnoticed  and  can  make  any  results  questionable.  Therefore,  we  decided  that  we  need  to  verify 
absolute  times  for  all  records  we  shall  use  in  the  future  before  we  attempt  anything  like  the 
separate  estimation  ot  compressional  and  deviatoric  source  energy. 
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FIGURE  41:  Surface  wave  imaging  of  two  nuclear  explosions  at  the  Kazakh  test 
site  (October  12,  1980  on  top  and  September  14,  1980  on  the  bottom)  by  the 
Capon  algorithm.  The  geographical  coordinates  of  the  centers  of  the  print  plot 
grids  are  50N*79E  and  the  events  are  plotted  relative  to  that.  Filled  circles  denote 
the  locations,  obtained  from  body  waves,  as  listed  by  unclassified  files  at  the 
Center  for  Seismic  Studies,  triangles  show  the  maxima  of  the  Capon’s  statistics. 
The  grids  are  50  km  on  each  side,  and  the  tops  are  facing  North. 
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Depth  imaging  of  the  1976  Kazakh  earthquake 


This  March  20,  1976  event  was  analyzed  in  some  detail  by  Pooley  et  al  (1983)  and  was 
identified  as  an  earthquake.  The  deconvolutions  at  three  UK  arrays  resulted  in  the  waveforms 
shown  in  Figure  42.  The  prominent  arrivals  include  P,  pP,  sP  and  an  early  arriving  Moho 
conversion  of  the  S  wave.  The  crustal  model  used  by  Pooley  et  al  (1983)  consisted  of  a  single 
layer  underlain  by  mantle.  Adapting  this  model  the  changes  in  the  depth  will  only  result  in  time 
shifts  of  the  three  idealized  impulsive  arrivals.  Using  the  Fourier  transforms  of  the  sequence  of 
these  four  impulses  spaced  appropriately  to  the  various  depths  as  /  and  those  of  the  three  data 
traces  for  d,  shown  in  Figure  43,  in  the  imaging  formula  we  can  see  how  the  frequency- 
compounded  F  statistic  between  0.5-4  Hz  changes  as  a  function  of  depth  (Figure  44).  We 
naturally  get  a  prominent  peak  at  the  expected  result  of  20.5  km,  having  assumed  the  same 
structure  and  phase  interpretation  as  Pooley  et  al  (1983).  This  is  an  obvious  example  with  clear 
depth  phases  that  could  be  interpreted  manually  as  well  and  it  merely  serves  as  an  illustration  how 
such  method  works.  We  do  not  claim  that  this  approach  would  work  so  neatly  all  the  time  for 
more  complicated  events  with  directionality  (Douglas  et  al,  1990).  We  must  point  out,  however, 
that  without  the  deconvolution  (or  equivalently,  beamed  phaseless  seismograms)  the  timings  and 
the  polarities  of  the  various  phases  would  be  much  less  obvious,  and  the  deconvolution  did  much 
to  clean  up  the  seismograms.  This  can  be  seen  in  an  earlier  repon  of  ours  where  we  presented 
samples  of  the  original  seismograms  for  this  event  (Der  et  al,  1987a). 

Summary  of  Section  III 

It  can  be  shown  that  time  domain  studies  of  long  period  body  waves  preferentially  fit  the 
low  frequency  components  of  signals  with  adverse  effect  on  the  resolution  of  the  spatial 
characteristics  of  the  seismic  sources  because  of  the  long  wavelengths  and  limited  resolution 
obtainable  from  long  period  signal  components.  Inversions  of  long  period  waveform  data  for 
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FIGURE  42:  Deconvolved  v^aveforms  of  the  Kazakh  earthquake  at  three  UK 
arrays.  Besides  the  surface  reflections  pP  and  sP  there  is  an  early  phase  which  is 
probably  a  Moho-converted  phase  (Pooley  et  al  1983). 
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FIGURE  44:  Comparisons  of  noise  adaptation  without  (bottom)  and  with  (top) 
the  site  response  removed.  Comparing  the  shapes  and  levels  of  spectra  it  can  be 
seen  that  the  overall  noise  is  signiHcantly  reduced  by  the  adaptive  beaming 
method  relatively  to  the  beams,  especially  below  2  Hz.  The  adaptive  beams 
quickly  adapt  to  the  noise,  and  after  a  few  oscillations  the  overall  noise  levels  are 
much  reduced  relative  to  the  beam  outputs.  EKA  data  was  processed. 


moderately  large  events  can  give  only  a  "centroid"  solution  without  any  reliable  details  in  the  fault 
motion.  Side-by-side  time  domain  comparisons  of  data  and  synthetic  seismograms  are  useful  aids 
for  verifying  the  results,  but  do  not  constitute  sufficient  proof  for  the  validity  of  the  source 
parameters  derived  by  such  comparisons  only. 

In  this  report,  we  present  alternative  general  ways  for  inverting  waveform  data.  The 
methods  are  a  general  class  of  pattern  matching  related  in  mathematical  form  to  F-K  analysis.  The 
methodology  has  the  advantage  that,  unlike  qualitative  comparison  of  data  and  synthetic 
waveforms,  it  provides  means  for  the  estimation  of  error  limits  and  spatial  resolution.  It  also 
appears  potentially  quit  suitable  for  automatization  of  the  waveform  inversion  process.  Preliminary 
applications  of  the  methodology  to  various  types  of  seismic  data,  long  period  and  short  period 
body  waves  and  long  period  surface  waves,  gave  interesting  and  promising  results.  A  potential 
application  for  such  algorithms  would  be  the  optimum  estimation  of  seismic  strain  release 
component  in  Ms-  Nevertheless,  much  needs  to  be  done  to  funher  develop  and  refine  such 
techniques,  especially  for  compound  sources. 
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4.0  TESTS  OF  VARIOUS  S/N  OPTIMINIZATION  ALGORITHMS 


General  Remarks 

In  applying  discrimination  and  yield  estimation  procedures  to  measurements  made  on 
various  teleseismic  and  regional  wave  arrivals  emitted  by  the  sources,  it  is  important  to  extend  the 
ranges  of  detectability  and  measurability.  Commonly,  large  events  will  excite  body  waves,  Lg, 
and  long  period  surface  waves  that  can  be  seen  teleseismically.  Smaller  events  will  not  have  long 
period  surface  waves  that  can  be  seen  above  the  noise  at  teleseismic  distances,  but  will  still  have 
detectable  P  waves.  Events  that  are  even  smaller  can  be  detected  only  by  regional  stations  and 
arrays.  With  decreasing  magnitude,  the  quantity  and  quality  of  measurements  used  as 
discriminants,  diagnostics  and  yield  estimators  will  also  decrease.  Therefore,  we  need  techniques 
that  extend  the  ranges  of  detectability  of  various  weak  arrivals  and  that  enable  us  to  extract 
broadband  information  from  the  data.  A  family  of  such  techniques  includes  the  various  optimum 
filtering  methods.  A  common  trait  of  many  of  these  is  that  they  were  designed  to  minimize 
stationary  noise  while  passing  ideal,  perfectly  identical  signal  waveforms  across  arrays  without 
distortion  (Capon,  1969;  Booker  and  Ong,  1971;  Shen,  1979;  Kedrov  and  Ovchinnikov,  1990). 

Past  attempts  of  optimum  filtering  were  characterized  by  performance  which  was  far  below 
that  predicted  from  idealized  theoretical  signal  and  noise  models  (Capon,  1969;  Booker  and  Ong, 
1971).  There  were  basically  two  reasons  for  this.  The  first  one  was  that  signal  waveforms  are  not 
exactly  identical  across  arrays,  this  led  to  degradation  (beam  loss)  in  the  performance  of  optimum 
filters  (Cox,  1973).  The  result  was  that  beam  loss  often  cancelled  the  gains  in  reducing  the  noise 
level  resulting  in  no  or  negligible  net  S/N  gain.  The  second  was  that  the  background  noise  is 
generally  not  stationary  and  the  noise  models  derived  from  large  data  blccks  were  not  realistic. 
Recognizing  the  nature  of  the  second  factor  several  workers  have  started  to  apply  adaptive  filtering 
methods  with  relatively  fast  adaptation  rate. 

Applying  such  techniques  to  small  arrays  where  the  signal  waveform  variability  was  not  a 
major  problem  significant  S/N  gains  were  indeed  achieved  (Widrow,  1966;  Shen,  1979).  He 
found  that  the  noise  reduction  exceeded  that  of  the  conventional  beamforming.  The  filters  were 
designed  on  the  basis  of  perfect  similarity  of  seismic  signals.  This  assumption  is  not  correct  for 
most  arrays,  as  we  have  pointed  out  above. 

In  this  report  we  are  investigating  the  effect  of  increasing  the  signal  similarity  by  factoring 
out  the  site  responses  prior  to  applying  the  technique  of  adaptive  filtering.  The  factoring  out  of  site 
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responses  will,  of  course,  change  the  noise  coherence  structure  also,  and  we  rely  on  the  adaptation 
of  the  filter  to  adjust  to  this  changed  structure.  After  removing  the  site  responses,  we  can  assume 
with  more  justification  that  the  signals  are  perfectly  coherent.  This  will  reduce  the  "beamloss" 
(Cox,  1973)  to  much  lower  levels. 

The  last  issue  is  that  of  the  computational  costs  associated  with  the  introduction  of  any  such 
relatively  elaborate  schemes.  In  the  past,  such  considerations  w'ere  very  important  because 
computers  were  expensive  and  slow  and  the  available  memory  or  disk  space  limited.  This  was  a 
major  reason,  besides  the  problems  discussed  above,  why  optimum  filtering  methods  fell  out  of 
favor  during  the  early  seventies.  With  the  advent  of  cheap  and  fast  computers  even  moderate  S/N 
gains  seem  to  be  worthwhile,  even  if  part  of  the  time  no  net  S/N  gains  result  relative  to  simpler 
processors.  It  can  be  shown  that  we  would  never  do  worse  than  a  simple  beam  processor  since 
optimum  processors  degrade  into  simple  beams  when  the  noise  is  uncorrelated.  Based  on  some 
limited  randomly  chosen  actual  signal  and  noise  situations  we  have  found,  however  that  net  S/N 
gains  (over  simple  beaming)  of  4-6  dB  are  quite  common.  Such  improvements  in  the  output  S/N 
seem  to  be  worthwhile. 

Processing  of  Teleseismic  P  Waves  Recorded  at  Large  Seismic  Arrays 

For  testing  our  approach,  we  have  used  recordings  of  a  set  of  Kazakh  (Shagan  River)  P 
waves  recorded  at  EKA  and  estimated  the  site  effects  by  applying  our  deconvolution  procedure  to 
the  seismograms  listed  in  Table  VIII.  Subsequently,  we  factored  out  the  site  effects  and  applied 
Shen’s  algorithm  to  the  result.  In  order  to  test  the  effect  of  the  algorithm  on  the  noise  alone,  we 
have  applied  the  noise-adaptive  filter  to  the  noise  preceding  the  signals. 

TABLE  VIII:  List  of  Kazakh  Events 
Recorded  at  AWRE  EKA  Array 


Event 

Origin  Time 

Lat  (N) 

Lon  (E) 

107C0/I  1 

A  >  /  A 

49.82 

78.14 

1978162 

02:56:57.6 

49.90 

78.80 

1977334 

04:06:57.4 

49.96 

78.89 

1976342 

04:56:57.4 

49.96 

78.85 
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We  performed  this  test  on  noise  samples  with  and  without  the  site  effect  correction. 
Figures  45  to  47  show  that  noise  adaptation  reduces  the  noise  amplitudes  relative  to  the 
conventional  beam  throughout  the  frequency  range  of  0-4  Hz,  but  the  noise  reduction  is  most 
effective  below  1  Hz.  The  removal  of  site  effects  did  not  impair  the  effectiveness  of  noise- 
adaptation.  The  filter  adapts  quickly  (as  evidenced  by  a  few  large  oscillations  in  the  output  after  the 
noise  amplitude  falls  off). 

Applying  the  process  to  windows  containing  signals  shows  that  the  noise  was  reduced  to  a 
level  below  visibility  in  this  high  S/N  sample  by  the  adaptive  process,  while  the  noise  is  clearly 
visible  on  the  beams  (top  traces  in  Figure  46).  Obviously,  this  example  has  no  practical  use,  since 
the  signal  is  much  larger  than  the  noise  in  the  case  shown.  It  must  be  pointed  out,  however,  that 
the  relative  noise  reduction  is  the  same  regardless  of  the  relative  signal  size  (the  site-equalized 
process  will  be  unaffected  by  the  presence  of  large  coherent  signals).  Thus,  a  weaker  signal  would 
have  stood  out  much  better  on  the  trace  processed  by  the  same  procedure  than  on  the  beam  and  the 
result  would  be  much  better  suited  for  any  further  analyses  of  the  signal. 

It  is  also  interesting  to  look  at  these  noise  reduction  results  in  the  frequency  domain.  In 
Figure  46  we  show  the  differences  in  the  output  on  signal  windows  when  the  adaptive  beam 
algorithm  is  applied  with  and  without  site  effect  correction.  The  figure  shows  that  without  the  site 
effect  correction  the  signal  output  is  considerably  less  with  adaptive  beaming  than  with 
conventional  beaming  (a).  This  is  not  surprising,  since  the  adaptive  beam  process  will  reduce  the 
signal  output  if  the  signals  are  not  matched.  After  site-equalization  is  performed,  this  signal  loss  is 
much  less  severe  (b).  Since  the  noise  reduction  of  adaptive  beaming  remains  more  effective,  even 
after  applying  site-equalization,  than  the  conventional  beam,  the  resulting  overall  S/N  ratio 
improvement  will  be  better  for  the  adaptive  beam  combined  with  site-equalization  (Figure  47). 
Thus,  this  combined  process  appears  to  be  superior  to  the  conventional  beaming  process.  The  gain 
is  achieved  at  a  considerable  expense  in  processing  complexity,  but  it  may  be  justified  in  special 
cases  when  we  want  to  learn  more  about  smaller  events  at  teleseismic  distances.  Such  procedures 
may  be  part  of  a  toolbox  for  more  detailed  investigations. 
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FIGURE  45:  Teleseismic  P  waves  recorded  at  EKA  processed  by  beaming  and 
adaptive  filtering  with  and  without  removal  of  site  frequency  responses. 
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FIGURE  46:  Spectra  representative  to  the  signal  beam  loss  (top)  and  the  S/N 
gain  obtained  for  two  events  by  using  the  four  modes  of  processing  described  in 
the  text. 


136 


5B0 


Sample  input  traces 


Conventionai  beam 


Adaptive  beam 


-.vvvA-v  VVy  Vt;  ly 
'\  A  A 


I  V  ' 


Conventional  beam 
with  site  corrections 


Adaptive  beam  with 
site  corrections 


I - 1 

3  sec 


FIGURE  47:  Optimum  processing  of  the  nuclear  explosion  Mast,  overlain  by 
added  amplified  actual  background  noise,  at  NORSAR.  The  best  performance, 
about  a  2:1  S/N  gain,  is  associated  with  site  equalization  followed  by  adaptive 
beaming  (bottom  trace).  It  must  be  pointed  out,  that  contrary  to  any  superficial 
appearances,  this  is  not  equivalent  to  frequency  filtering  since  this  process  works 
by  spatial  filtering  only  and  does  not  change  the  spectrum  or  waveforms  of  any 
signals. 
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Another  example  of  this  kind  of  signal  processing  is  that  of  the  P  waves  from  the  NTS 
nuclear  explosion  Mast  recorded  at  NORSAR.  To  simulate  noisy  conditions,  we  have  amplified 
background  noise  preceding  the  event  and  added  it  to  the  respective  traces  of  the  array  recording. 
The  site  response  functions  were  obtained  by  factoring  a  suite  of  original,  high  S/N  recordings  of 
NTS  nuclear  explosions  at  NORSAR.  Thus,  we  are  simulating  the  case  of  enhancing  a  weak  event 
by  utilizing  pre-stored  high-quality  site  response  functions  in  conjunction  with  adaptive  filtering. 
In  this  simulation,  we  have  utilized  12  sensors  of  NORSAR  spread  out  across  the  reduced  array. 
Figure  48  shows  some  of  the  original,  unprocessed  traces  on  the  top  and  the  results  of  four  kinds 
of  processing.  It  is  clear  that  the  combination  of  the  site-equalization  and  the  adaptive  noise 
filtering  gave  the  best  result. 

Processing  of  Western  Norway  Earthquakes  Recorded  at  NORESS 

Similarly  to  teleseismic  signals,  we  have  done  a  limited  amount  of  work  testing  the 
suitability  of  adaptive  beaming  and  site-equalization  methods  for  enhancing  regional  Pn  arrivals. 
Pn  arrivals  are  generally  emergent  at  regional  distances  and  are  much  smaller  in  amplitude  than  Lg 
or  Sn.  Enhancing  Pn  arrivals  may  be  quite  useful  in  discrimination  studies.  The  methodology 
followed  is  identical  to  that  described  for  teleseismic  P  elsewhere  in  this  repon.  The  method  was 
tested  using  a  set  of  Pn  arrivals  from  western  Norway  earthquakes.  A  list  of  these  events  is  given 
in  Table  IX.  The  site  effects  were  estimated  from  applying  the  factorization  procedure  using  the 
multi-channel  deconvolution  program  (Der  et  al,  1987).  Applying  the  adaptive  beaming  program 
to  the  noise  alone  (Figure  48)  shows  that  the  site  response  removal  somewhat  degrades  the  relative 
performance  improvement  due  to  the  adaptive  beaming  process.  On  the  other  hand,  the  decrease  in 
signal  loss  through  the  process  appears  to  compensate  approximately  for  this  loss.  The  net  result 
is  that  the  best  strategy  seems  to  be  to  apply  adaptive  (not  conventional)  beaming  either  with  or 
without  site  corrections.  We  must  add  that  this  test  was  done  before  we  realized  the  extreme 
sensitivity  of  site  factors  to  azimuth,  and  the  western  Norway  earthquakes  may  be  too  scattered  in 
azimuth  to  get  stable  and  effective  site  correction  factors,  in  any  case,  the  adaptive  beaming 
process  appears  to  be  better  than  conventional  beaming  for  regional  signals. 
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1  98521  7 
Norway  data 
Noise  beams 


1 98521 7 
Norway  data 
Noise  beams 
Site  resp  removed 


1985217 
Norway  data 
Signal  beams 


1985217 
Norway  data 
Signal  beams 
Site  resp  removed 


FIGURE  48:  Noise  reduction  gain  and  beam  loss  for  western  Norway 

earthquakes  as  processed  by  beaming  and  adaptive  filtering  with  or  without 
removing  the  site  transfer  functions. 
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TABLE  IX;  List  of  Western  Norway  Events  Recorded  at  NORESS 


Event 

Origin  Time 

Lat  (N) 

Lon  (E) 

1985217 

17:42:58.7 

59.3 

6.59 

1985218 

17:50:07.9 

59.3 

6.95 

1985290 

10:00:00.4 

59.3 

6.95 

-  CONCLUSIONS 

Work  during  this  contract  revealed  the  great  potential  in  the  applications  of  various 
statistical  time-series  analysis  methods  for  solving  seismological  problems  relevant  to  nuclear 
monitoring.  In  our  study  of  the  coherence  structure  regional  seismic  waveforms,  extending  the 
concept  of  spectral  factorability  of  teleseismic  body  wave  data  (Filson  and  Frasier,  1972),  we  have 
demonstrated  that  spectra  of  regional  phases  can  also  be  decomposed  into  source  and  recording  site 
factors  provided  that  the  signals  originated  from  a  limited  source  region.  Groups  of  events  located 
close  to  each  other  can  be  identified  by  the  fact  that  their  waveforms  can  be  reconstructed  from  their 
spectral  source  and  site  factors. 

In  addition,  it  was  also  found  that  events  within  such  factorable  groups  could  be  further 
subdivided  according  to  relative  inter-event  coherences  between  them.  The  events  which  show 
high  inter-event  coherence  are  probably  both  close  to  each  other  and  have  very  similar  source 
mechanisms,  albeit  different  source  time  functions.  The  opposite  must  be  true  for  events  that  do 
not  show  high  coherence.  Cross-event  equalization  filtering  between  coherent  events  using  short 
time  domain  filters  resulted  in  increases  in  waveform  similarity,  and  no  such  increase  seemed 
possible  between  events  pairs  with  low  inter-event  coherence. 

Site  transfer  functions  between  sensors  located  at  the  extreme  ends  of  NORESS  appear  to 
be  complex,  describable  only  with  transfer  functions  with  impulse  responses  longer  than  a  second. 
Assuming  the  same  degree  of  crustal  heterogeneity  for  a  source  region  this  may  imply  changes  in 
waveforms  similar  in  nature  over  small  displacements  of  source  location. 
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Analyses  of  the  coherence  properties  of  suites  of  regional  waveforms  for  Pn  and  Lg  can  be 
used  for: 

a)  Identifying  very  closely  spaced  events  with  similar  source  mechanisms. 

b)  Identifying  closely  spaced  events  with  different  source  mechanisms  (possibly 
attempts  of  evasion). 

c)  Identifying  larger  spatial  separations  in  events,  not  detectable  by  standard 
location  methods. 

In  the  second  part  of  this  report,  we  re-examined  some  aspects  of  estimating  secondary 
arrivals  (pP,  spall,  etc.)  in  teleseismic  P  waves.  A  popular  model  (P+pP  model)  often  used  for 
interpreting  P  waves  from  nuclear  explosions  and  deriving  pP  parameters  consist  of  a  site  function 
convolved  with  a  P  and  pP  pulse  sequence  and  some  explosion  source  pulse  (Lay,  1985;  Murphy, 
1989).  Frequency  domain  analyses  of  data  variance  show  that  this  model  cannot  explain  a 
significant  portion  of  the  energy  (30-40%)  in  the  P  waves  from  Pahute  Mesa  events  for  which  such 
analysis  methods  were  previously  applied.  Spectral  factorization,  which  does  not  assume  the 
P+pP  model,  on  the  other  hand,  accounts  for  90-95%  of  the  energy  in  both  Pahute  Mesa  and 
Kazakh  events.  Therefore,  the  complex  source  waveforms  often  seen  in  multi-channel 
deconvolved  P  waves  from  nuclear  explosions  at  both  Pahute  Mesa  and  the  Kazakh  test  sites  are 
required  for  explaining  the  data  and  reflect  the  actual  complexity  of  P  waves  radiated  from  nuclear 
explosions.  It  is  shown  that  magnitude-yield  relationships  appropriate  to  the  generic  pP 
characteristics  derived  from  multi-channel  deconvolutions  of  nuclear  explosions  (and  appropriate 
t*)  are  similar  to  those  directly  derived  by  regression  analyses  of  mb  data  and  published  yields  of 
nuclear  explosions  for  Pahute  and  Kazakh  test  sites  by  Jih  (1990). 

In  the  third  section  of  this  report,  we  have  derived  some  new  approaches  for  imaging 
seismic  sources  and  seismic  source  inversion.  The  methods  are  based  finding  the  maxima  of 
statistics  that  are  closely  related  to  those  used  in  F-K  analyses  as  functions  of  source  parameters. 
The  method  is  quite  general,  it  uses  Green's  functions  in  conjunction  with  waveforms  of  various 
types  of  arrivals.  Such  can  be  used  for  inverting  both  body  and  surface  waves  for  source 
parameters  in  a  semi-automatic  fashion  without  subjective  visual  or  RMS  comparisons  of  synthetic 
and  data  waveforms.  These  algorithms  also  provide  a  means  for  evaluation  of  the  performance  of 
time  domain  direct  modelling  and  source  inversion  methods  as  currently  practiced.  Our  analyses 
indicated  that  the  time  domain  fits  using  long  period  data  are  generally  not  precise  enough  to 
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resolve  the  source  depth  to  a  few  km  and  other  spatial  details  with  the  accuracies  often  claimed,  but 
can  only  give  "centroid"  solutions.  Other  simulations  indicated  that  the  same  approaches  can  create 
reasonably  resolved  images  of  faults  and  can  be  used  for  the  separate  estimation  of  compressional 
and  strain  release  components  of  energies  from  explosions. 

Our  source  imaging-inversion  methodology  was  applied  to  two  deconvolved  Kazakh 
nuclear  explosions  and  the  1966  El  Golfo,  the  1988  Armenian  and  1976  Kazakh  earthquakes.  The 
method  identified  both  explosions  as  shallow,  with  depths  below  5  km,  and  the  depth  estimate  for 
the  El  Golfo  and  Kazakh  earthquakes  event  turned  out  to  be  14  km  and  20.5  km.  These  are 
reasonable  results  indicating  that  such  methods  could  be  applied  for  seismic  discrimination. 
Imaging  the  Armenian  earthquake  gave  a  linear  feature  close  in  orientation  to  the  actual  fault 
inferred  as  the  source  of  this  event.  Imaging  Kazakh  explosion  sources  using  surface  waves  was 
hampered  by  demonstrably  erroneous  absolute  timing  information  in  the  CSS  data  we  retrieved. 

In  the  last  section  of  this  report,  we  explored  applications  of  combined  adaptive  filtering 
and  site-effect  compensation  processing  schemes  to  teleseismic  (UK  array)  and  regional 
(NORESS)  data  indicated  that  such  methods  can  significantly  enhance  the  signal- to-noise  ratios  (by 
4-6  dB)  for  P  (Pn)  arrivals  should  such  computer-intensive  analyses  prove  to  be  desirable  for  an 
event. 
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APPENDIX  A 


Source  imaging:  A  statistical  theory 

As  we  stated  in  the  introduction  of  this  report,  one  of  the  main  goals  of  analysis  of  seismic 
data  is  to  derive  some  properties  of  the  source.  In  nuclear  monitoring  we  desire  to  discriminate 
between  events  of  various  types,  such  as  earthquakes  and  explosions  either  through  the  direct 
discrimination  of  source  mechanisms  or  indirectly  through  determination  of  the  source  depth.  It 
would  be  desirable  to  do  this  by  some  visual  identification  of  sources  after  some  preprocessing.  In 
conventional  processing  of  P  wave  data  such  identification  can  be  performed  through  lining  up  the 
short  period  P  waves  and  observing  the  pP  "moveout",  i.e.,  the  changes  in  pP  time  with  the  values 
of  slowness.  Alternatively,  one  may  apply  some  criteria  using  amplitude  measurements  in  some 
time  windows  covering  possible  times  of  depth  phases  to  ascertain  that  an  event  belong  to  either 
the  earthquake  or  explosion  category  (Pearce  1980).  Although  such  methods  are  clearly  useful, 
they  are  somewhat  tedious,  and  some  future  automation  of  such  approaches  is  clearly  desirable.  In 
a  broader  sense  going  beyond  the  concept  of  point  sources,  detailed  mapping  of  the  sources  of 
seismic  radiation  is  also  very  important  in  studies  of  the  mechanisms  of  earthquakes.  In  order  to 
understand  the  faulting  processes  we  need  to  know  how  the  fault  dislocation  developed  in  time  and 
space. 


The  discussion  below  establishes  a  proposed  definition  of  the  "source  image",  i.e.,  the 
estimated  power  of  seismic  sources  as  functions  of  space  coordinates  in  the  source  region  and 
derives  confidence  intervals  for  this  image.  We  consider,  as  usual,  the  observed  time  vector  y  (t) 
observed  over  N  array  elements  or  stations.  The  discrete  Fourier  transform  of  this  can  be  written 
as 


y  =f  s  +  n 


(1) 


A-1 


where  s  =  s(x,  is  the  Fourier  transform  of  the  theoretical  signal  at  frequency  vand  at  source 
position  vector  ;c,  usually  specified  as  a  combination  of  location  and  depth.  The  vector/ =  f(x,  v) 
is  the  transform  of  the  function  that  combines  with  s  to  produce  the  data  y.  We  can  think  of  /  as 
some  combination  of  Green's  functions.  The  vector/is  assumed  to  be  known  whereas  the  signal  s 
is  deterministic  and  unknown.  The  noise  vector  n  is  assumed  to  be  spatially  white  (uncorrelated 
among  stations)  with  power  spectral  matrix 

£(« /i*)  =  />„(v)/  (2)j 

at  frequency  v,  where  Pn(  v'  J  is  the  noise  power  spectrum  of  each  channel  and  1  is  the  identity 
matrix. 


Information  carried  by  the  observed  vector  y.  One  can  simply  plot  the  power  in  the 
generalized  transform/"  y  as  the  beampower 


B{x,v)  = 


An  alternative  measure  with  higher  resolution  is  Capon's  C,  defined  as 


C  (jc,  V )  = 


1/ 


•(y 


y  + 


c^I 


)■' 


/ 


-1 


(3) 


(4) 


which  can  also  be  written  in  the  form 
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which  shows  its  relation  to  the  generalized  beam.  It  is  convenient  to  define  the  residual  power  as 


v )  =  [y  -  B  (x,  V )  (6) 

so  that 

I  B  ( X  V)  \ 

C(x,v)=-^1  + - .  (7) 

[/■I  (c^  +  /?(x,  V)) 

/ 

the  constant  c  is  an  arbitrary  number  chosen  to  make  the  matrix  (yy*  +  I)  non-singular. 

A  second  alternative  estimator  is  Shumway’s  F,  defined  as  the  likelihood  ratio  test  statistic 
resulting  from  testing  the  model  (1).  This  leads  to 

which  is  closely  related  to  the  other  two. 

Capon's  C  and  Shumway's  F  have  similar  resolving  capabilities  since  they  are  both 
functions  of  B  and  R.  They  are  both  better  than  the  generalized  beampower  from  this  point  of 
view  since  equation  (3)  does  not  incorporate  the  residual  noise  power  R.  We  shall  show  below 
that  the  F  is  superior  statistically  to  the  other  two,  both  for  detecting  the  predefined  signal  and  for 
providing  an  estimator  for  the  signal  image. 


If  the  noise  is  Gaussian,  it  is  easy  to  show  that 

P^iv) 

(9) 

t 

where  x^2  denotes  the  chi-squared  distribution  with  noncentrality  parameter 

2  2l  SU.vf 

Pn(V) 

(10) 

and  ~  is  a  notation  for  "approximately  distributed  as".  Note  that  5^  is  twice  the 

"signal-to-noise" 

ratio 

^  Pn(y) 

(11) 

It  is  this  signal-to-noise  ratio  that  we  propose  to  define  as  the  image.  It  is  also  easy  to  show  that 

2  R  {x,v)  2  /  Q  s 

P„iv) 

(12) 

and  is  independent  of  B(x,  v).  It  follows  that 

• 

r  ,s\  xl^sSn 

^2,2N-2^°  9 

xl,.2li2N-2) 

(13) 

SO  that 

(14) 
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has  the  non-central  F  distribution  with  the  non-centrality  parameter  5^  given  by  equation  (11). 
Hence  the  distribution  of  F  depends  only  on  the  signal-to-noise  ratio  ^(x,  v)  for  a  given  number 
of  sensors  N. 


We  may  examine  this  further  by  noting  that  for  an  F  with  nj  and  «2  degrees  of  freedom 


"l."2 


”2 


^2  ■  2 


.  (rtj  +  5  ) 


(15) 


so  that 


E{F(x,v))  =  (1  +  )) 


(16) 


and  the  mean  value  is  essentially  the  signal-to-noise  rauo. 


‘2  “/V'‘2 


can  be  used  lO  compute  the  variance  with 


var  iF)  =  E{F^)  -E\f) 


(18) 


and  nj  =  2,  np  =  2  N  -  2.  The  variance  depends  on  the  signal-to-noise  ratio  that  we  are  trying  to 
estimate  in  a  rather  complicated  way  and  the  normal  approximation  will  not  be  very  good  for 
obtaining  an  approximate  confidence  interval. 
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For  any  jc  and  v,  however,  we  can  obtain  an  exact  confidence  interval  for  E(F(x,v)),  8i^  or  in 
the  following  way.  We  know  from  (16)  that 


^  (x,v) 


(N-2) 

(N-l) 


F{x,v)-\ 


(19) 


is  an  unbiased  estimator  of  the  signal  to  noise  ratio  (11)  which  we  are  defining  as  the  image.  Let 
F(x,v)  be  some  observed  sample  value  of  equation  (14).  Let  and  6^2  i^o  values  of 

such  that 


Pr  {  2^5^  )<F  {x,v)  }  =  a/2 


(20) 


and 


{  ^2.2N-2(^2^  )<F{x,v)  )  =  a/2  (21) 

The  resulting  1-a  confidence  interval  for  ^^will  be  (5^j,  S^2)-  The  interval  for  this 
signal-to-noise  ratio  is  (1/2  i,  1/2  3^2)-  This  interval  can  also  be  connected  into  a  (1-a) 
confidence  interval  for  E(F(x,v))  using  equation  (16). 

The  above  procedure  can  be  used  to  get  a  (1-a)  percent  confidence  interval  for  each  x  and 
V.  Since  the  intervals  are  dependent,  for  n  of  them  one  can  only  assign  an  overall  confidence  of  1- 
na.  Therefore,  in  order  to  assign  small  confidence  to  the  upper  and  lower  surfaces  of,  say  .90, 
with  100  points  ion  the  surface,  each  separate  point  should  have  a=.001  or  99.9%  confidence. 

One  may  also  investigate  various  pooling  procedures  if  the  model  is  repeated  at  M  points, 
say 
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(22) 


for  m-l,2,M,  where  fm  =  f  (x,  Vm),  Sm  =  s  (x,  Vm)  aiid  the  spectrum  of  the  noise  is 


(22) 


so  that  it  is  relatively  constant  over  the  frequencies.  Then,  for 


B  {X,  , 


W‘ 


(24) 


(25) 


where 


g2_2|5(x,v^)f 
p  tv) 


It  follows  that  the  pooled  value 


(26) 


M 


M 


^'Lb^x.v^) 

p  ^  y  j  ^2JH  ( ( 


m=l 


m=\ 


(27) 


Also,  for  the  residual  power 


(28) 
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it  follows  that 


M 


2£ 


m=l 


Pni^) 


■  ~ XlM {N-\) 


Hence,  we  have  the  result  that 


F{x,v^)  =  {N-\)^ 


B  (X,  vj 


is  distributed  as  a  non-central  F  with  2M  and  2  M  (N-1 )  degrees  of  freedom 
1978)  and  non-centrality  parameter 


M 


2  X  v„  )f 


S  = 


m=l 


^(V') 


One  might  also  consider  pooling  the  separate  F  values 


leading  to 


F(x,v^)  =  iN-l) 


B  (x,  vj 
R  (X,  v^) 


M 


m=\ 
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(29) 


(30) 

(Hogg  and  Craig 


(31) 


(32) 


(33) 


An  approximate  interval  might  be  based  on  normality  and  homogeneity  (not  likely  to  be  true 
because  of  (17))  of  the  means  and  variances  over  Vm*  Then  a  very  rough  and  ready  approximation 
would  be 


F2  ix,  v^)±z 


(34) 


where 


M 

m=l 


and  Z(x/2  is  the  upper  a/2  point  on  the  normal  distribution. 
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APPENDIX  B 


An  alternative  method  for  imaging 

Although  the  estimation  of  source  spatio-temporal  configurations  can  be  done  effectively 
using  the  adaptation  of  Shumway’s  F-K  method  computing  of  confidence  limits  becomes 
somewhat  crude  if  we  assume  normality  and  if  we  have  only  a  few  frequencies  that  contribute  to 
the  image.  A  normal  distribution  in  this  case  is  a  very  poor  approximation  to  the  superposition  of 
noncentral  F  distributions.  To  derive  an  alternative  approach  we  again  start  out  with  the  model  in 
the  frequency  domain,  assuming  that  all  quantities  are  complex  and  functions  of  frequency 

y  =f  s  +  n  (1) 

where  the  nx  1  probe  (Green's  function)  vector / relates  the  scalar  signal  s  to  the  n x  ? data  vector 
with  an  additive  vector  of  noises  n  assumed  to  be  complex  Gaussian  with  zero  mean  and  a 
complex  covariance  matrix 


E{nn*)  =  PJ  (2) 

where  /  is  the  n  x  n  identity  matrix.  We  assume  that  the  scalar  signal  s  is  univariate  complex 
Gaussian  with  power  density 


E{s*s)  =  P^  (3) 

which  we  are  trying  to  estimate  as  a  function  of  /’s  to  form  an  image.  Assuming  that  it  is 
stochastic  rather  then  fixed  makes  the  derivation  of  statistics  easier,  we  need  to  work  with  only 
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central  chi-squared  variables.  The  following  arguments  are  essentially  adaptations  of  those  in 
chapter  4.4.5  of  Shumway  (1988). 


The  key  is  to  organize  the  test  statistics  that  we  have  looked  at  before  as  a  spectral  analysis 
of  variance  problem.  Denote  the  numerator  B  of  our  F  statistics  in  the  previous  section  as 

B  ^  f  (4) 

// 

and  interpret  it  as  the  power  due  to  the  signal  since  it  is  essentially  the  extra  power  (over  tr(5) 
where  S  is  the  spectral  matrix  of  the  data)  accounted  for  in  fitting  the  signal  model.  Denote  the 
denominator  by 


R  =  tr(S)-B  (5) 

which  is  the  noise  or  residual  error  power. 

This  can  be  arranged  in  an  analysis  of  power  table  as  follows 

Source  Power  Degrees  of  freedom  E(power) 


Signal 


fSf* 

ff* 


2  ff*Ps+Pn 


Noise 


f  f  * 


2  (N  -  1) 


(N-1)P„ 
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Total 


tr(S) 


2N 


Where  N  is  the  number  of  stations  applied.  We  verify  the  expectation  by  noting  that  the 
beam  power  B  is  approximately  distributed  as 


X\2) 

2 


Then,  since 


(6) 


we  get  the  result.  Note  also  that 


var 


x\m) 


=  2m 


(7) 


so  that 


var 


+  P. 

s  n 


(8) 


R  is  approximately  distributed  as 


P,x\2N  -2) 
2 


so  that  E(R)  =  Pfi  and 


var  iR)  =  iN  -1)  P^ 
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Now  consider  an  unbiased  estimator  for  the  image  Ps,  defined  as 


B - - 

(N  -1) 

// 


(9) 


It  is  clear  that  ^  ^  and  that 


var  (P^) 


(10) 


It  is  clear  that  both  the  estimator  and  the  estimator  for  its  variance  can  be  computed  for  the 
analysis  of  power  table  but  we  need  to  define  the  confidence  intervals. 

To  obtain  the  confidence  interval  we  use  an  old  stratagem  due  originally  to  Sattertwaite 
(Ref?)  and  also  used  by  Tukey  (Ref?)  to  approximate  the  distribution  of  smoother  spectral 
estimates  with  windows.  We  assume  that  the  quantity  which  is  the  weighted  difference  of  chi- 
square  variable  can  also  be  approximated  with  a  chi-squared  variable.  We  shall  need  the  degrees  of 
freedom  for  this  chi-squared  distribution.  This  can  be  estimated  by  taking  the  ratio  of  twice  the 
expectation  to  the  variance. 


2[E\x\m)^ 
var  [x\m  )] 


(11) 


In  our  case  this  becomes 
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m  = 


(12) 


var  (P^ ) 


where  ^has  been  computed  earlier.  The  confidence  interval  can  be  computed  using  the  fact 

that 


has  a  distribution. 


Naturally,  it  is  necessary  to  replace  all  quantities  that  are  unknown  in  these  expressions 
with  their  estimators  from  the  analysis  of  power  table  above. 


To  derive  confidence  intervals  for  signal  power  estimated  from  pooling  L  independent 
values  for  various  frequencies,  the  degrees  of  freedom  m'  for  the  pooled  statistics,  also  assumed  to 
be  distributed  as  chi-squared,  we  get  the  degrees  of  freedom  from  the  expression. 


m 


i=l 


(13) 
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APPENDIX  C 


Implementation  of  multi-channel  deconvolution  in  the  X  windows  environment 

In  order  to  facilitate  the  use  and  the  testing  of  the  multi-channel  deconvolution  method  we 
have  written  a  version  that  makes  use  of  the  X-window  environment.  The  program  is  menu  driven 
with  the  following  options  for  viewing  in  the  Multi-channel  Deconvolution  Display  (MDD): 

1)  Source  estimate  waveforms  and  resolution  kernels. 

2)  Site  transfer  function  (time  domain)  displays. 

3)  Source  spectral  estimates. 

4)  Site  response  spectra. 

5)  Original  vs.  reconstructed  waveforms  with  associated  correlation  coefficients. 

Examples  of  these  various  kinds  of  displays  are  shown  in  figures  Cl  tt'  C4.  For  further 
explanations,  we  refer  to  the  figure  captions.  When  MDD  is  activated  it  creates  a  window  and 
waits  for  user  input.  All  input  is  entered  by  using  a  mouse.  The  right  mouse  button  produces  a 
menu  of  the  options.  The  option  is  chosen  by  moving  the  mouse  so  that  the  desired  option  is 
highlighted  and  the  left  mouse  button  is  pressed. 

Once  the  option  is  chosen  a  window  displaying  the  desired  data  appears.  If  all  the  data  will 
not  fit  the  window  the  first  few  waveforms  will  be  displayed  and  the  rest  can  be  viewed  by 
scrolling.  A  scrollbar  appears  at  the  left  edge  of  each  window.  It  can  be  used  to  scroll  the  data 
within  that  window  by  standard  X 1 1  scrollbar  manipulations.  Each  new  window  will  be  created 
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on  top  of  those  currently  displayed.  All  windows  can  be  resized,  moved  and  overlayed  using 
standard  X  facilities. 
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Multichannel  Deconvolut 
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